Intrepid
example_16.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38 // Denis Ridzal (dridzal@sandia.gov), or
39 // Kara Peterson (kjpeter@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
81 // Intrepid includes
85 #include "Intrepid_CellTools.hpp"
90 #include "Intrepid_Utils.hpp"
91 
92 // Epetra includes
93 #include "Epetra_Time.h"
94 #include "Epetra_Map.h"
95 #include "Epetra_FEVector.h"
96 #include "Epetra_SerialComm.h"
97 
98 // Teuchos includes
99 #include "Teuchos_oblackholestream.hpp"
100 #include "Teuchos_RCP.hpp"
101 #include "Teuchos_Time.hpp"
102 
103 // Shards includes
104 #include "Shards_CellTopology.hpp"
105 
106 // EpetraExt includes
107 #include "EpetraExt_MultiVectorOut.h"
108 
109 using namespace std;
110 using namespace Intrepid;
111 using Teuchos::rcp;
112 
113 int main(int argc, char *argv[]) {
114 
115  //Check number of arguments
116  if (argc < 4) {
117  std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
118  std::cout <<"Usage:\n\n";
119  std::cout <<" ./Intrepid_example_Drivers_Example_14.exe deg NX NY NZ verbose\n\n";
120  std::cout <<" where \n";
121  std::cout <<" int deg - polynomial degree to be used (assumed >= 1) \n";
122  std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n";
123  std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n";
124  std::cout <<" int NZ - num intervals in y direction (assumed box domain, 0,1) \n";
125  std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
126  exit(1);
127  }
128 
129  // This little trick lets us print to std::cout only if
130  // a (dummy) command-line argument is provided.
131  int iprint = argc - 1;
132  Teuchos::RCP<std::ostream> outStream;
133  Teuchos::oblackholestream bhs; // outputs nothing
134  if (iprint > 2)
135  outStream = Teuchos::rcp(&std::cout, false);
136  else
137  outStream = Teuchos::rcp(&bhs, false);
138 
139  // Save the format state of the original std::cout.
140  Teuchos::oblackholestream oldFormatState;
141  oldFormatState.copyfmt(std::cout);
142 
143  *outStream \
144  << "===============================================================================\n" \
145  << "| |\n" \
146  << "| Example: Apply Stiffness Matrix for |\n" \
147  << "| Poisson Equation on Hexahedral Mesh |\n" \
148  << "| |\n" \
149  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
150  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
151  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
152  << "| |\n" \
153  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
154  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
155  << "| |\n" \
156  << "===============================================================================\n";
157 
158 
159  // ************************************ GET INPUTS **************************************
160 
161  int deg = atoi(argv[1]); // polynomial degree to use
162  int NX = atoi(argv[2]); // num intervals in x direction (assumed box domain, 0,1)
163  int NY = atoi(argv[3]); // num intervals in y direction (assumed box domain, 0,1)
164  int NZ = atoi(argv[4]); // num intervals in y direction (assumed box domain, 0,1)
165 
166 
167  // *********************************** CELL TOPOLOGY **********************************
168 
169  // Get cell topology for base hexahedron
170  typedef shards::CellTopology CellTopology;
171  CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() );
172 
173  // Get dimensions
174  int numNodesPerElem = hex_8.getNodeCount();
175  int spaceDim = hex_8.getDimension();
176 
177  // *********************************** GENERATE MESH ************************************
178 
179  *outStream << "Generating mesh ... \n\n";
180 
181  *outStream << " NX" << " NY" << " NZ\n";
182  *outStream << std::setw(5) << NX <<
183  std::setw(5) << NY << std::setw(5) << NZ << "\n\n";
184 
185  // Print mesh information
186  int numElems = NX*NY*NZ;
187  int numNodes = (NX+1)*(NY+1)*(NZ+1);
188  *outStream << " Number of Elements: " << numElems << " \n";
189  *outStream << " Number of Nodes: " << numNodes << " \n\n";
190 
191  // Cube
192  double leftX = 0.0, rightX = 1.0;
193  double leftY = 0.0, rightY = 1.0;
194  double leftZ = 0.0, rightZ = 1.0;
195 
196  // Mesh spacing
197  double hx = (rightX-leftX)/((double)NX);
198  double hy = (rightY-leftY)/((double)NY);
199  double hz = (rightZ-leftZ)/((double)NZ);
200 
201  // Get nodal coordinates
202  FieldContainer<double> nodeCoord(numNodes, spaceDim);
203  FieldContainer<int> nodeOnBoundary(numNodes);
204  int inode = 0;
205  for (int k=0; k<NZ+1; k++)
206  {
207  for (int j=0; j<NY+1; j++)
208  {
209  for (int i=0; i<NX+1; i++)
210  {
211  nodeCoord(inode,0) = leftX + (double)i*hx;
212  nodeCoord(inode,1) = leftY + (double)j*hy;
213  nodeCoord(inode,2) = leftZ + (double)k*hz;
214  if (k==0 || k==NZ || j==0 || i==0 || j==NY || i==NX)
215  {
216  nodeOnBoundary(inode)=1;
217  }
218  else
219  {
220  nodeOnBoundary(inode)=0;
221  }
222  inode++;
223  }
224  }
225  }
226 //#define DUMP_DATA
227 #ifdef DUMP_DATA
228  // Print nodal coords
229  ofstream fcoordout("coords.dat");
230  for (int i=0; i<numNodes; i++) {
231  fcoordout << nodeCoord(i,0) <<" ";
232  fcoordout << nodeCoord(i,1) <<" ";
233  fcoordout << nodeCoord(i,2) <<"\n";
234  }
235  fcoordout.close();
236 #endif
237 
238 
239 
240  // ********************************* CUBATURE AND BASIS ***********************************
241  *outStream << "Getting cubature and basis ... \n\n";
242 
243  // Get numerical integration points and weights
244  // I only need this on the line since I'm doing tensor products
245  Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > glcub
246  = Teuchos::rcp(new CubaturePolylib<double,FieldContainer<double>,FieldContainer<double> >(2*deg-1,PL_GAUSS_LOBATTO) );
247 
248  const int numCubPoints1D = glcub->getNumPoints();
249 
250  FieldContainer<double> cubPoints1D(numCubPoints1D, 1);
251  FieldContainer<double> cubWeights1D(numCubPoints1D);
252 
253  glcub->getCubature(cubPoints1D,cubWeights1D);
254 
255  std::vector<Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > >
256  cub_to_tensor;
257  cub_to_tensor.push_back( glcub );
258  cub_to_tensor.push_back( glcub );
259  cub_to_tensor.push_back( glcub );
260 
261  Array<RCP<FieldContainer<double> > > wts_by_dim(3);
262  wts_by_dim[0] = rcp( &cubWeights1D , false ); wts_by_dim[1] = wts_by_dim[0]; wts_by_dim[2] = wts_by_dim[1];
263 
265 
266  Basis_HGRAD_HEX_Cn_FEM<double, FieldContainer<double> > hexBasis( deg , POINTTYPE_SPECTRAL );
267 
268  Array< Array< RCP< Basis< double , FieldContainer<double> > > > > &bases = hexBasis.getBases();
269 
270  // get the bases tabulated at the quadrature points, dimension-by-dimension
271 
272  Array< RCP< FieldContainer<double> > > basisVals( 3 );
273  FieldContainer<double> bvals1D( bases[0][0]->getCardinality() , numCubPoints1D );
274  bases[0][0]->getValues( bvals1D , cubPoints1D , OPERATOR_VALUE );
275  basisVals[0] = rcp( &bvals1D , false ); basisVals[1] = basisVals[0]; basisVals[2] = basisVals[0];
276 
277  Array< RCP< FieldContainer<double> > > basisDVals( 3 );
278  FieldContainer<double> bdvals1D( bases[0][0]->getCardinality() , numCubPoints1D , 1);
279  bases[0][0]->getValues( bdvals1D , cubPoints1D , OPERATOR_D1 );
280  basisDVals[0] = rcp( &bdvals1D , false ); basisDVals[1] = basisDVals[0]; basisDVals[2] = basisDVals[0];
281 
282 
283  const int numCubPoints = cubhex.getNumPoints();
284  FieldContainer<double> cubPoints3D( numCubPoints , 3 );
285  FieldContainer<double> cubWeights3D( numCubPoints );
286  cubhex.getCubature( cubPoints3D , cubWeights3D );
287 
288 
289  FieldContainer<int> elemToNode(numElems, numNodesPerElem);
290  int ielem = 0;
291  for (int k=0; k<NZ; k++)
292  {
293  for (int j=0; j<NY; j++)
294  {
295  for (int i=0; i<NX; i++)
296  {
297  elemToNode(ielem,0) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
298  elemToNode(ielem,1) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
299  elemToNode(ielem,2) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
300  elemToNode(ielem,3) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
301  elemToNode(ielem,4) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
302  elemToNode(ielem,5) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
303  elemToNode(ielem,6) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
304  elemToNode(ielem,7) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
305  ielem++;
306  }
307  }
308  }
309 #ifdef DUMP_DATA
310  // Output connectivity
311  ofstream fe2nout("elem2node.dat");
312  for (int k=0;k<NZ;k++)
313  {
314  for (int j=0; j<NY; j++)
315  {
316  for (int i=0; i<NX; i++)
317  {
318  int ielem = i + j * NX + k * NY * NY;
319  for (int m=0; m<numNodesPerElem; m++)
320  {
321  fe2nout << elemToNode(ielem,m) <<" ";
322  }
323  fe2nout <<"\n";
324  }
325  }
326  }
327  fe2nout.close();
328 #endif
329 
330 
331  // ********************************* 3-D LOCAL-TO-GLOBAL MAPPING *******************************
332  FieldContainer<int> ltgMapping(numElems,hexBasis.getCardinality());
333  const int numDOF = (NX*deg+1)*(NY*deg+1)*(NZ*deg+1);
334  ielem=0;
335  for (int k=0;k<NZ;k++)
336  {
337  for (int j=0;j<NY;j++)
338  {
339  for (int i=0;i<NX;i++)
340  {
341  const int start = k * ( NY * deg + 1 ) * ( NX * deg + 1 ) + j * ( NX * deg + 1 ) + i * deg;
342  // loop over local dof on this cell
343  int local_dof_cur=0;
344  for (int kloc=0;kloc<=deg;kloc++)
345  {
346  for (int jloc=0;jloc<=deg;jloc++)
347  {
348  for (int iloc=0;iloc<=deg;iloc++)
349  {
350  ltgMapping(ielem,local_dof_cur) = start
351  + kloc * ( NX * deg + 1 ) * ( NY * deg + 1 )
352  + jloc * ( NX * deg + 1 )
353  + iloc;
354  local_dof_cur++;
355  }
356  }
357  }
358  ielem++;
359  }
360  }
361  }
362 #ifdef DUMP_DATA
363  // Output ltg mapping
364  ielem = 0;
365  ofstream ltgout("ltg.dat");
366  for (int k=0;k<NZ;k++)
367  {
368  for (int j=0; j<NY; j++)
369  {
370  for (int i=0; i<NX; i++)
371  {
372  int ielem = i + j * NX + k * NX * NY;
373  for (int m=0; m<hexBasis.getCardinality(); m++)
374  {
375  ltgout << ltgMapping(ielem,m) <<" ";
376  }
377  ltgout <<"\n";
378  }
379  }
380  }
381  ltgout.close();
382 #endif
383 
384  // ********** DECLARE GLOBAL OBJECTS *************
385  Epetra_SerialComm Comm;
386  Epetra_Map globalMapG(numDOF, 0, Comm);
387  Epetra_FEVector u(globalMapG); u.Random();
388  Epetra_FEVector Ku(globalMapG);
389 
390  // ************* For Jacobians **********************
391  FieldContainer<double> cellVertices(numElems,numNodesPerElem,spaceDim);
392  FieldContainer<double> cellJacobian(numElems,numCubPoints,spaceDim,spaceDim);
393  FieldContainer<double> cellJacobInv(numElems,numCubPoints,spaceDim,spaceDim);
394  FieldContainer<double> cellJacobDet(numElems,numCubPoints);
395 
396 
397  // get vertices of cells (for computing Jacobians)
398  for (int i=0;i<numElems;i++)
399  {
400  for (int j=0;j<numNodesPerElem;j++)
401  {
402  const int nodeCur = elemToNode(i,j);
403  for (int k=0;k<spaceDim;k++)
404  {
405  cellVertices(i,j,k) = nodeCoord(nodeCur,k);
406  }
407  }
408  }
409 
410  // jacobian evaluation
411  CellTools<double>::setJacobian(cellJacobian,cubPoints3D,cellVertices,hex_8);
412  CellTools<double>::setJacobianInv(cellJacobInv, cellJacobian );
413  CellTools<double>::setJacobianDet(cellJacobDet, cellJacobian );
414 
415 
416  // ************* MATRIX-FREE APPLICATION
417  FieldContainer<double> uScattered(numElems,1,hexBasis.getCardinality());
418  FieldContainer<double> KuScattered(numElems,1,hexBasis.getCardinality());
419  FieldContainer<double> gradU(numElems,1,hexBasis.getCardinality(),3);
420 
421  u.GlobalAssemble();
422 
423 
424 
425  Ku.PutScalar(0.0);
426  Ku.GlobalAssemble();
427 
428  double *uVals = u[0];
429  double *KuVals = Ku[0];
430 
431  Teuchos::Time full_timer( "Time to apply operator matrix-free:" );
432  Teuchos::Time scatter_timer( "Time to scatter dof:" );
433  Teuchos::Time elementwise_timer( "Time to do elementwise computation:" );
434  Teuchos::Time grad_timer( "Time to compute gradients:" );
435  Teuchos::Time pointwise_timer( "Time to do pointwise transformations:" );
436  Teuchos::Time moment_timer( "Time to compute moments:" );
437  Teuchos::Time gather_timer( "Time to gather dof:" );
438 
439  full_timer.start();
440 
441  scatter_timer.start();
442  for (int k=0; k<numElems; k++)
443  {
444  for (int i=0;i<hexBasis.getCardinality();i++)
445  {
446  uScattered(k,0,i) = uVals[ltgMapping(k,i)];
447  }
448  }
449  scatter_timer.stop();
450 
451  elementwise_timer.start();
452 
453  grad_timer.start();
454  Intrepid::TensorProductSpaceTools::evaluateGradient<double>( gradU , uScattered ,basisVals , basisDVals );
455  grad_timer.stop();
456  pointwise_timer.start();
457  Intrepid::FunctionSpaceToolsInPlace::HGRADtransformGRAD<double>( gradU , cellJacobian );
458  Intrepid::FunctionSpaceToolsInPlace::HGRADtransformGRADDual<double>( gradU , cellJacobian );
459  Intrepid::FunctionSpaceToolsInPlace::multiplyMeasure<double>( gradU , cellJacobDet );
460  pointwise_timer.stop();
461  moment_timer.start();
462  Intrepid::TensorProductSpaceTools::momentsGrad<double>( KuScattered , gradU , basisVals , basisDVals , wts_by_dim );
463  moment_timer.stop();
464  elementwise_timer.stop();
465  gather_timer.start();
466  for (int k=0;k<numElems;k++)
467  {
468  for (int i=0;i<hexBasis.getCardinality();i++)
469  {
470  KuVals[ltgMapping(k,i)] += KuScattered(k,0,i);
471  }
472  }
473  gather_timer.stop();
474  full_timer.stop();
475 
476  *outStream << full_timer.name() << " " << full_timer.totalElapsedTime() << " sec\n";
477  *outStream << "\t" << scatter_timer.name() << " " << scatter_timer.totalElapsedTime() << " sec\n";
478  *outStream << "\t" << elementwise_timer.name() << " " << elementwise_timer.totalElapsedTime() << " sec\n";
479  *outStream << "\t\t" << grad_timer.name() << " " << grad_timer.totalElapsedTime() << " sec\n";
480  *outStream << "\t\t" << pointwise_timer.name() << " " << pointwise_timer.totalElapsedTime() << " sec\n";
481  *outStream << "\t\t" << moment_timer.name() << " " << moment_timer.totalElapsedTime() << " sec\n";
482  *outStream << "\t" << gather_timer.name() << " " << gather_timer.totalElapsedTime() << " sec\n";
483 
484 
485  *outStream << "End Result: TEST PASSED\n";
486 
487  // reset format state of std::cout
488  std::cout.copyfmt(oldFormatState);
489 
490  return 0;
491 }
492 
Header file for the Intrepid::FunctionSpaceToolsInPlace class.
Header file for the Intrepid::CellTools class.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::CubatureTensor class.
Header file for the Intrepid::HGRAD_HEX_Cn_FEM class.
Intrepid utilities.
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin, Aeronautics, Imperial College London) within Intrepid.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell...
Header file for the Intrepid::FunctionSpaceTools class.
Header file for the Intrepid::CubaturePolylib class.
Defines tensor-product cubature (integration) rules in Intrepid.
Header file for the Intrepid::TensorProductSpaceTools class.
A stateless class for operations on cell data. Provides methods for: