Intrepid
example_13.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 
82 // Intrepid includes
85 #include "Intrepid_CellTools.hpp"
86 //#include "Intrepid_ArrayTools.hpp"
88 //#include "Intrepid_RealSpaceTools.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_BLAS.hpp"
102 //#include "Teuchos_BLAS_types.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 
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_10.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: Build 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  // Element to Node map
240  // We'll keep it around, but this is only the DOFMap if you are in the lowest order case.
241  FieldContainer<int> elemToNode(numElems, numNodesPerElem);
242  int ielem=0;
243  for (int k=0; k<NZ; k++)
244  {
245  for (int j=0; j<NY; j++)
246  {
247  for (int i=0; i<NX; i++)
248  {
249  elemToNode(ielem,0) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
250  elemToNode(ielem,1) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
251  elemToNode(ielem,2) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
252  elemToNode(ielem,3) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
253  elemToNode(ielem,4) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
254  elemToNode(ielem,5) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
255  elemToNode(ielem,6) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
256  elemToNode(ielem,7) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
257  ielem++;
258  }
259  }
260  }
261 #ifdef DUMP_DATA
262  // Output connectivity
263  ofstream fe2nout("elem2node.dat");
264  for (int k=0;k<NZ;k++)
265  {
266  for (int j=0; j<NY; j++)
267  {
268  for (int i=0; i<NX; i++)
269  {
270  ielem = i + j * NX + k * NY * NY;
271  for (int m=0; m<numNodesPerElem; m++)
272  {
273  fe2nout << elemToNode(ielem,m) <<" ";
274  }
275  fe2nout <<"\n";
276  }
277  }
278  }
279  fe2nout.close();
280 #endif
281 
282  // ************************************ CUBATURE **************************************
283  *outStream << "Getting cubature ... \n\n";
284 
285  // Get numerical integration points and weights
286  DefaultCubatureFactory<double> cubFactory;
287  int cubDegree = 2*deg;
288  Teuchos::RCP<Cubature<double> > quadCub = cubFactory.create(hex_8, cubDegree);
289 
290  int cubDim = quadCub->getDimension();
291  int numCubPoints = quadCub->getNumPoints();
292 
293  FieldContainer<double> cubPoints(numCubPoints, cubDim);
294  FieldContainer<double> cubWeights(numCubPoints);
295 
296  quadCub->getCubature(cubPoints, cubWeights);
297 
298  // ************************************** BASIS ***************************************
299 
300  *outStream << "Getting basis ... \n\n";
301 
302  // Define basis
303  Basis_HGRAD_HEX_Cn_FEM<double, FieldContainer<double> > quadHGradBasis(deg,POINTTYPE_SPECTRAL);
304  int numFieldsG = quadHGradBasis.getCardinality();
305  FieldContainer<double> quadGVals(numFieldsG, numCubPoints);
306  FieldContainer<double> quadGrads(numFieldsG, numCubPoints, spaceDim);
307 
308  // Evaluate basis values and gradients at cubature points
309  quadHGradBasis.getValues(quadGVals, cubPoints, OPERATOR_VALUE);
310  quadHGradBasis.getValues(quadGrads, cubPoints, OPERATOR_GRAD);
311 
312  // create the local-global mapping
313  FieldContainer<int> ltgMapping(numElems,numFieldsG);
314  const int numDOF = (NX*deg+1)*(NY*deg+1)*(NZ*deg+1);
315  ielem=0;
316  for (int k=0;k<NZ;k++)
317  {
318  for (int j=0;j<NY;j++)
319  {
320  for (int i=0;i<NX;i++)
321  {
322  const int start = k * ( NY * deg + 1 ) * ( NX * deg + 1 ) + j * ( NX * deg + 1 ) + i * deg;
323  // loop over local dof on this cell
324  int local_dof_cur=0;
325  for (int kloc=0;kloc<=deg;kloc++)
326  {
327  for (int jloc=0;jloc<=deg;jloc++)
328  {
329  for (int iloc=0;iloc<=deg;iloc++)
330  {
331  ltgMapping(ielem,local_dof_cur) = start
332  + kloc * ( NX * deg + 1 ) * ( NY * deg + 1 )
333  + jloc * ( NX * deg + 1 )
334  + iloc;
335  local_dof_cur++;
336  }
337  }
338  }
339  ielem++;
340  }
341  }
342  }
343 #ifdef DUMP_DATA
344  // Output ltg mapping
345  ofstream ltgout("ltg.dat");
346  for (int k=0;k<NZ;k++)
347  {
348  for (int j=0; j<NY; j++)
349  {
350  for (int i=0; i<NX; i++)
351  {
352  ielem = i + j * NX + k * NX * NY;
353  for (int m=0; m<numFieldsG; m++)
354  {
355  ltgout << ltgMapping(ielem,m) <<" ";
356  }
357  ltgout <<"\n";
358  }
359  }
360  }
361  ltgout.close();
362 #endif
363 
364  // ********** DECLARE GLOBAL OBJECTS *************
365  Epetra_SerialComm Comm;
366  Epetra_Map globalMapG(numDOF, 0, Comm);
367  Epetra_FEVector u(globalMapG); u.Random();
368  Epetra_FEVector Ku(globalMapG);
369 
370 
371 
372  // ********** CONSTRUCT AND INSERT LOCAL STIFFNESS MATRICES ***********
373  *outStream << "Building reference stiffness matrix...\n\n";
375  typedef FunctionSpaceTools fst;
376 
377  // jacobian information
378  FieldContainer<double> refCellNodes(1,numNodesPerElem,spaceDim);
379  FieldContainer<double> cellJacobian(1,numCubPoints,spaceDim,spaceDim);
380  FieldContainer<double> cellJacobInv(1,numCubPoints,spaceDim,spaceDim);
381  FieldContainer<double> cellJacobDet(1,numCubPoints);
382 
383  // element stiffness matrices and supporting storage space
384  FieldContainer<double> localStiffMatrix(1, numFieldsG, numFieldsG);
385  FieldContainer<double> transformedBasisGradients(1,numFieldsG,numCubPoints,spaceDim);
386  FieldContainer<double> weightedTransformedBasisGradients(1,numFieldsG,numCubPoints,spaceDim);
387  FieldContainer<double> weightedMeasure(1, numCubPoints);
388 
389  Epetra_Time localConstructTimer( Comm );
390  refCellNodes(0,0,0) = 0.0; refCellNodes(0,0,1) = 0.0; refCellNodes(0,0,2) = 0.0;
391  refCellNodes(0,1,0) = hx; refCellNodes(0,1,1) = 0.0; refCellNodes(0,1,2) = 0.0;
392  refCellNodes(0,2,0) = hx; refCellNodes(0,2,1) = hy; refCellNodes(0,2,2) = 0.0;
393  refCellNodes(0,3,0) = 0.0; refCellNodes(0,3,1) = hy; refCellNodes(0,3,2) = 0.0;
394  refCellNodes(0,4,0) = 0.0; refCellNodes(0,4,1) = 0.0; refCellNodes(0,4,2) = hz;
395  refCellNodes(0,5,0) = hx; refCellNodes(0,5,1) = 0.0; refCellNodes(0,5,2) = hz;
396  refCellNodes(0,6,0) = hx; refCellNodes(0,6,1) = hy; refCellNodes(0,6,2) = hz;
397  refCellNodes(0,7,0) = 0.0; refCellNodes(0,7,1) = hy; refCellNodes(0,7,2) = hz;
398 
399  // jacobian evaluation
400  CellTools::setJacobian(cellJacobian,cubPoints,refCellNodes,hex_8);
401  CellTools::setJacobianInv(cellJacobInv, cellJacobian );
402  CellTools::setJacobianDet(cellJacobDet, cellJacobian );
403 
404  // transform reference element gradients to each cell
405  fst::HGRADtransformGRAD<double>(transformedBasisGradients, cellJacobInv, quadGrads);
406 
407  // compute weighted measure
408  fst::computeCellMeasure<double>(weightedMeasure, cellJacobDet, cubWeights);
409 
410  // multiply values with weighted measure
411  fst::multiplyMeasure<double>(weightedTransformedBasisGradients,
412  weightedMeasure, transformedBasisGradients);
413 
414  // integrate to compute element stiffness matrix
415  fst::integrate<double>(localStiffMatrix,
416  transformedBasisGradients, weightedTransformedBasisGradients , COMP_BLAS);
417 
418  const double localConstructTime = localConstructTimer.ElapsedTime();
419 
420  // ************* MATRIX-FREE APPLICATION
421  FieldContainer<double> uScattered(numElems,numFieldsG);
422  FieldContainer<double> KuScattered(numElems,numFieldsG);
423 
424  u.GlobalAssemble();
425 
426  Epetra_Time multTimer(Comm);
427  Teuchos::BLAS<int,double> blas;
428  Ku.PutScalar(0.0);
429  Ku.GlobalAssemble();
430 
431  double *uVals = u[0];
432  double *KuVals = Ku[0];
433 
434  Epetra_Time scatterTimer(Comm);
435  std::cout << "Scattering\n";
436  // Scatter
437  for (int k=0; k<numElems; k++)
438  {
439  for (int i=0;i<numFieldsG;i++)
440  {
441  uScattered(k,i) = uVals[ltgMapping(k,i)];
442  }
443  }
444  const double scatterTime = scatterTimer.ElapsedTime();
445 
446  Epetra_Time blasTimer(Comm);
447  blas.GEMM(Teuchos::NO_TRANS , Teuchos::NO_TRANS ,
448  numFieldsG , numElems, numFieldsG ,
449  1.0 ,
450  &localStiffMatrix(0,0,0) ,
451  numFieldsG ,
452  &uScattered(0,0) ,
453  numFieldsG ,
454  0.0 ,
455  &KuScattered(0,0) ,
456  numFieldsG );
457  const double blasTime = blasTimer.ElapsedTime();
458 
459  Epetra_Time gatherTimer(Comm);
460  // Gather
461  for (int k=0;k<numElems;k++)
462  {
463  for (int i=0;i<numFieldsG;i++)
464  {
465  KuVals[ltgMapping(k,i)] += KuScattered(k,i);
466  }
467  }
468 
469  const double gatherTime = gatherTimer.ElapsedTime();
470 
471 
472  *outStream << "Time to build local matrix (including Jacobian computation): "<< localConstructTime << "\n";
473  *outStream << "Time to scatter " << scatterTime << "\n";
474  *outStream << "Time for local application " << blasTime << "\n";
475  *outStream << "Time to gather " << gatherTime << "\n";
476  *outStream << "Total matrix-free time " << scatterTime + blasTime + gatherTime << "\n";
477 
478 
479  *outStream << "End Result: TEST PASSED\n";
480 
481  // reset format state of std::cout
482  std::cout.copyfmt(oldFormatState);
483 
484  return 0;
485 }
486 
Header file for the Intrepid::CellTools class.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_HEX_Cn_FEM class.
Defines expert-level interfaces for the evaluation of functions and operators in physical space (supp...
Intrepid utilities.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell...
Header file for the Intrepid::FunctionSpaceTools class.
A factory class that generates specific instances of cubatures.
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > &degree)
Factory method.
A stateless class for operations on cell data. Provides methods for: