Intrepid
example_17.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 #include "Teuchos_SerialDenseMatrix.hpp"
103 
104 // Shards includes
105 #include "Shards_CellTopology.hpp"
106 
107 // EpetraExt includes
108 #include "EpetraExt_MultiVectorOut.h"
109 
110 using namespace std;
111 using namespace Intrepid;
112 using Teuchos::rcp;
113 
114 int main(int argc, char *argv[]) {
115 
116  //Check number of arguments
117  if (argc < 4) {
118  std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
119  std::cout <<"Usage:\n\n";
120  std::cout <<" ./Intrepid_example_Drivers_Example_14.exe deg NX NY NZ verbose\n\n";
121  std::cout <<" where \n";
122  std::cout <<" int deg - polynomial degree to be used (assumed >= 1) \n";
123  std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n";
124  std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n";
125  std::cout <<" int NZ - num intervals in y direction (assumed box domain, 0,1) \n";
126  std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
127  exit(1);
128  }
129 
130  // This little trick lets us print to std::cout only if
131  // a (dummy) command-line argument is provided.
132  int iprint = argc - 1;
133  Teuchos::RCP<std::ostream> outStream;
134  Teuchos::oblackholestream bhs; // outputs nothing
135  if (iprint > 2)
136  outStream = Teuchos::rcp(&std::cout, false);
137  else
138  outStream = Teuchos::rcp(&bhs, false);
139 
140  // Save the format state of the original std::cout.
141  Teuchos::oblackholestream oldFormatState;
142  oldFormatState.copyfmt(std::cout);
143 
144  *outStream \
145  << "===============================================================================\n" \
146  << "| |\n" \
147  << "| Example: Apply Stiffness Matrix for |\n" \
148  << "| Poisson Operator on Hexahedral Mesh with BLAS |\n" \
149  << "| |\n" \
150  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
151  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
152  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
153  << "| |\n" \
154  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
155  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
156  << "| |\n" \
157  << "===============================================================================\n";
158 
159 
160  // ************************************ GET INPUTS **************************************
161 
162  int deg = atoi(argv[1]); // polynomial degree to use
163  int NX = atoi(argv[2]); // num intervals in x direction (assumed box domain, 0,1)
164  int NY = atoi(argv[3]); // num intervals in y direction (assumed box domain, 0,1)
165  int NZ = atoi(argv[4]); // num intervals in y direction (assumed box domain, 0,1)
166 
167 
168  // *********************************** CELL TOPOLOGY **********************************
169 
170  // Get cell topology for base hexahedron
171  typedef shards::CellTopology CellTopology;
172  CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() );
173 
174  // Get dimensions
175  int numNodesPerElem = hex_8.getNodeCount();
176  int spaceDim = hex_8.getDimension();
177 
178  // *********************************** GENERATE MESH ************************************
179 
180  *outStream << "Generating mesh ... \n\n";
181 
182  *outStream << " NX" << " NY" << " NZ\n";
183  *outStream << std::setw(5) << NX <<
184  std::setw(5) << NY << std::setw(5) << NZ << "\n\n";
185 
186  // Print mesh information
187  int numElems = NX*NY*NZ;
188  int numNodes = (NX+1)*(NY+1)*(NZ+1);
189  *outStream << " Number of Elements: " << numElems << " \n";
190  *outStream << " Number of Nodes: " << numNodes << " \n\n";
191 
192  // Cube
193  double leftX = 0.0, rightX = 1.0;
194  double leftY = 0.0, rightY = 1.0;
195  double leftZ = 0.0, rightZ = 1.0;
196 
197  // Mesh spacing
198  double hx = (rightX-leftX)/((double)NX);
199  double hy = (rightY-leftY)/((double)NY);
200  double hz = (rightZ-leftZ)/((double)NZ);
201 
202  // Get nodal coordinates
203  FieldContainer<double> nodeCoord(numNodes, spaceDim);
204  FieldContainer<int> nodeOnBoundary(numNodes);
205  int inode = 0;
206  for (int k=0; k<NZ+1; k++)
207  {
208  for (int j=0; j<NY+1; j++)
209  {
210  for (int i=0; i<NX+1; i++)
211  {
212  nodeCoord(inode,0) = leftX + (double)i*hx;
213  nodeCoord(inode,1) = leftY + (double)j*hy;
214  nodeCoord(inode,2) = leftZ + (double)k*hz;
215  if (k==0 || k==NZ || j==0 || i==0 || j==NY || i==NX)
216  {
217  nodeOnBoundary(inode)=1;
218  }
219  else
220  {
221  nodeOnBoundary(inode)=0;
222  }
223  inode++;
224  }
225  }
226  }
227 //#define DUMP_DATA
228 #ifdef DUMP_DATA
229  // Print nodal coords
230  ofstream fcoordout("coords.dat");
231  for (int i=0; i<numNodes; i++) {
232  fcoordout << nodeCoord(i,0) <<" ";
233  fcoordout << nodeCoord(i,1) <<" ";
234  fcoordout << nodeCoord(i,2) <<"\n";
235  }
236  fcoordout.close();
237 #endif
238 
239 
240 
241  // ********************************* CUBATURE AND BASIS ***********************************
242  *outStream << "Getting cubature and basis ... \n\n";
243 
244  // Get numerical integration points and weights
245  // I only need this on the line since I'm doing tensor products
246  Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > glcub
247  = Teuchos::rcp(new CubaturePolylib<double,FieldContainer<double>,FieldContainer<double> >(2*deg-1,PL_GAUSS_LOBATTO) );
248 
249  const int numCubPoints1D = glcub->getNumPoints();
250 
251  FieldContainer<double> cubPoints1D(numCubPoints1D, 1);
252  FieldContainer<double> cubWeights1D(numCubPoints1D);
253 
254  glcub->getCubature(cubPoints1D,cubWeights1D);
255 
256  std::vector<Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > >
257  cub_to_tensor;
258  cub_to_tensor.push_back( glcub );
259  cub_to_tensor.push_back( glcub );
260  cub_to_tensor.push_back( glcub );
261 
263 
264  Basis_HGRAD_HEX_Cn_FEM<double, FieldContainer<double> > hexBasis( deg , POINTTYPE_SPECTRAL );
265  const int numFields = hexBasis.getCardinality();
266 
267  // get the bases tabulated at the quadrature points, dimension-by-dimension
268  const int numCubPoints = cubhex.getNumPoints();
269  FieldContainer<double> cubPoints3D( numCubPoints , spaceDim );
270  FieldContainer<double> cubWeights3D( numCubPoints );
271  FieldContainer<double> basisGrads( numFields , numCubPoints , spaceDim );
272  cubhex.getCubature( cubPoints3D , cubWeights3D );
273  hexBasis.getValues( basisGrads , cubPoints3D , OPERATOR_GRAD );
274 
275 
276  FieldContainer<int> elemToNode(numElems, numNodesPerElem);
277  int ielem = 0;
278  for (int k=0; k<NZ; k++)
279  {
280  for (int j=0; j<NY; j++)
281  {
282  for (int i=0; i<NX; i++)
283  {
284  elemToNode(ielem,0) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
285  elemToNode(ielem,1) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
286  elemToNode(ielem,2) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
287  elemToNode(ielem,3) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
288  elemToNode(ielem,4) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
289  elemToNode(ielem,5) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
290  elemToNode(ielem,6) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
291  elemToNode(ielem,7) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
292  ielem++;
293  }
294  }
295  }
296 #ifdef DUMP_DATA
297  // Output connectivity
298  ofstream fe2nout("elem2node.dat");
299  for (int k=0;k<NZ;k++)
300  {
301  for (int j=0; j<NY; j++)
302  {
303  for (int i=0; i<NX; i++)
304  {
305  ielem = i + j * NX + k * NY * NY;
306  for (int m=0; m<numNodesPerElem; m++)
307  {
308  fe2nout << elemToNode(ielem,m) <<" ";
309  }
310  fe2nout <<"\n";
311  }
312  }
313  }
314  fe2nout.close();
315 #endif
316 
317 
318  // ********************************* 3-D LOCAL-TO-GLOBAL MAPPING *******************************
319  FieldContainer<int> ltgMapping(numElems,numFields);
320  const int numDOF = (NX*deg+1)*(NY*deg+1)*(NZ*deg+1);
321  ielem=0;
322  for (int k=0;k<NZ;k++)
323  {
324  for (int j=0;j<NY;j++)
325  {
326  for (int i=0;i<NX;i++)
327  {
328  const int start = k * ( NY * deg + 1 ) * ( NX * deg + 1 ) + j * ( NX * deg + 1 ) + i * deg;
329  // loop over local dof on this cell
330  int local_dof_cur=0;
331  for (int kloc=0;kloc<=deg;kloc++)
332  {
333  for (int jloc=0;jloc<=deg;jloc++)
334  {
335  for (int iloc=0;iloc<=deg;iloc++)
336  {
337  ltgMapping(ielem,local_dof_cur) = start
338  + kloc * ( NX * deg + 1 ) * ( NY * deg + 1 )
339  + jloc * ( NX * deg + 1 )
340  + iloc;
341  local_dof_cur++;
342  }
343  }
344  }
345  ielem++;
346  }
347  }
348  }
349 #ifdef DUMP_DATA
350  // Output ltg mapping
351  ielem = 0;
352  ofstream ltgout("ltg.dat");
353  for (int k=0;k<NZ;k++)
354  {
355  for (int j=0; j<NY; j++)
356  {
357  for (int i=0; i<NX; i++)
358  {
359  ielem = i + j * NX + k * NX * NY;
360  for (int m=0; m<numFields; m++)
361  {
362  ltgout << ltgMapping(ielem,m) <<" ";
363  }
364  ltgout <<"\n";
365  }
366  }
367  }
368  ltgout.close();
369 #endif
370 
371  // ********** DECLARE GLOBAL OBJECTS *************
372  Epetra_SerialComm Comm;
373  Epetra_Map globalMapG(numDOF, 0, Comm);
374  Epetra_FEVector u(globalMapG); u.Random();
375  Epetra_FEVector Ku(globalMapG);
376 
377  // ************* For Jacobians **********************
378  FieldContainer<double> cellVertices(numElems,numNodesPerElem,spaceDim);
379  FieldContainer<double> cellJacobian(numElems,numCubPoints,spaceDim,spaceDim);
380  FieldContainer<double> cellJacobInv(numElems,numCubPoints,spaceDim,spaceDim);
381  FieldContainer<double> cellJacobDet(numElems,numCubPoints);
382 
383 
384  // get vertices of cells (for computing Jacobians)
385  for (int i=0;i<numElems;i++)
386  {
387  for (int j=0;j<numNodesPerElem;j++)
388  {
389  const int nodeCur = elemToNode(i,j);
390  for (int k=0;k<spaceDim;k++)
391  {
392  cellVertices(i,j,k) = nodeCoord(nodeCur,k);
393  }
394  }
395  }
396 
397  // jacobian evaluation
398  CellTools<double>::setJacobian(cellJacobian,cubPoints3D,cellVertices,hex_8);
399  CellTools<double>::setJacobianInv(cellJacobInv, cellJacobian );
400  CellTools<double>::setJacobianDet(cellJacobDet, cellJacobian );
401 
402 
403  // ************* MATRIX-FREE APPLICATION
404  FieldContainer<double> uScattered(numElems,1,numFields);
405  FieldContainer<double> KuScattered(numElems,1,numFields);
406  FieldContainer<double> gradU(numElems,1,numFields,3);
407 
408  // get views in SDM format for BLAS purposes
409  Teuchos::SerialDenseMatrix<int,double> refGradsAsMat( Teuchos::View ,
410  &basisGrads[0] ,
411  numCubPoints * spaceDim ,
412  numCubPoints * spaceDim ,
413  numFields );
414  Teuchos::SerialDenseMatrix<int,double> uScatAsMat( Teuchos::View ,
415  &uScattered[0] ,
416  numFields ,
417  numFields ,
418  numElems );
419  Teuchos::SerialDenseMatrix<int,double> GraduScatAsMat( Teuchos::View ,
420  &gradU[0] ,
421  numCubPoints * spaceDim ,
422  numCubPoints * spaceDim ,
423  numElems );
424  Teuchos::SerialDenseMatrix<int,double> KuScatAsMat( Teuchos::View ,
425  &KuScattered[0] ,
426  numFields ,
427  numFields ,
428  numElems );
429 
430 
431  u.GlobalAssemble();
432 
433 
434 
435  Ku.PutScalar(0.0);
436  Ku.GlobalAssemble();
437 
438  double *uVals = u[0];
439  double *KuVals = Ku[0];
440 
441  Teuchos::Time full_timer( "Time to apply operator matrix-free:" );
442  Teuchos::Time scatter_timer( "Time to scatter dof:" );
443  Teuchos::Time elementwise_timer( "Time to do elementwise computation:" );
444  Teuchos::Time grad_timer( "Time to compute gradients:" );
445  Teuchos::Time pointwise_timer( "Time to do pointwise transformations:" );
446  Teuchos::Time moment_timer( "Time to compute moments:" );
447  Teuchos::Time gather_timer( "Time to gather dof:" );
448 
449  full_timer.start();
450 
451  scatter_timer.start();
452  for (int k=0; k<numElems; k++)
453  {
454  for (int i=0;i<numFields;i++)
455  {
456  uScattered(k,0,i) = uVals[ltgMapping(k,i)];
457  }
458  }
459  scatter_timer.stop();
460 
461  elementwise_timer.start();
462 
463  grad_timer.start();
464  // reference grad per cell with BLAS
465  GraduScatAsMat.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS ,
466  1.0 , refGradsAsMat , uScatAsMat , 0.0 );
467 
468 
469  grad_timer.stop();
470  pointwise_timer.start();
471  // pointwise transformations of gradients
472  Intrepid::FunctionSpaceToolsInPlace::HGRADtransformGRAD<double>( gradU , cellJacobian );
473  Intrepid::FunctionSpaceToolsInPlace::HGRADtransformGRADDual<double>( gradU , cellJacobian );
474  Intrepid::FunctionSpaceToolsInPlace::multiplyMeasure<double>( gradU , cellJacobDet );
475  pointwise_timer.stop();
476  moment_timer.start();
477 
478  KuScatAsMat.multiply( Teuchos::TRANS , Teuchos::NO_TRANS , 1.0 , refGradsAsMat , GraduScatAsMat , 0.0 );
479 
480  moment_timer.stop();
481  elementwise_timer.stop();
482  gather_timer.start();
483  for (int k=0;k<numElems;k++)
484  {
485  for (int i=0;i<numFields;i++)
486  {
487  KuVals[ltgMapping(k,i)] += KuScattered(k,0,i);
488  }
489  }
490  gather_timer.stop();
491  full_timer.stop();
492 
493  *outStream << full_timer.name() << " " << full_timer.totalElapsedTime() << " sec\n";
494  *outStream << "\t" << scatter_timer.name() << " " << scatter_timer.totalElapsedTime() << " sec\n";
495  *outStream << "\t" << elementwise_timer.name() << " " << elementwise_timer.totalElapsedTime() << " sec\n";
496  *outStream << "\t\t" << grad_timer.name() << " " << grad_timer.totalElapsedTime() << " sec\n";
497  *outStream << "\t\t" << pointwise_timer.name() << " " << pointwise_timer.totalElapsedTime() << " sec\n";
498  *outStream << "\t\t" << moment_timer.name() << " " << moment_timer.totalElapsedTime() << " sec\n";
499  *outStream << "\t" << gather_timer.name() << " " << gather_timer.totalElapsedTime() << " sec\n";
500 
501 
502  *outStream << "End Result: TEST PASSED\n";
503 
504  // reset format state of std::cout
505  std::cout.copyfmt(oldFormatState);
506 
507  return 0;
508 }
509 
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: