Intrepid
example_03.cpp
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 
80 // Intrepid includes
83 #include "Intrepid_CellTools.hpp"
84 #include "Intrepid_ArrayTools.hpp"
85 #include "Intrepid_HGRAD_HEX_C1_FEM.hpp"
88 #include "Intrepid_Utils.hpp"
89 
90 // Epetra includes
91 #include "Epetra_Time.h"
92 #include "Epetra_Map.h"
93 #include "Epetra_FECrsMatrix.h"
94 #include "Epetra_FEVector.h"
95 #include "Epetra_SerialComm.h"
96 
97 // Teuchos includes
98 #include "Teuchos_oblackholestream.hpp"
99 #include "Teuchos_RCP.hpp"
100 #include "Teuchos_BLAS.hpp"
101 
102 // Shards includes
103 #include "Shards_CellTopology.hpp"
104 
105 // EpetraExt includes
106 #include "EpetraExt_RowMatrixOut.h"
107 #include "EpetraExt_MultiVectorOut.h"
108 
109 using namespace std;
110 using namespace Intrepid;
111 
112 // Functions to evaluate exact solution and derivatives
113 double evalu(double & x, double & y, double & z);
114 int evalGradu(double & x, double & y, double & z, double & gradu1, double & gradu2, double & gradu3);
115 double evalDivGradu(double & x, double & y, double & z);
116 
117 int main(int argc, char *argv[]) {
118 
119  //Check number of arguments
120  if (argc < 4) {
121  std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
122  std::cout <<"Usage:\n\n";
123  std::cout <<" ./Intrepid_example_Drivers_Example_03.exe NX NY NZ verbose\n\n";
124  std::cout <<" where \n";
125  std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n";
126  std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n";
127  std::cout <<" int NZ - num intervals in z direction (assumed box domain, 0,1) \n";
128  std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
129  exit(1);
130  }
131 
132  // This little trick lets us print to std::cout only if
133  // a (dummy) command-line argument is provided.
134  int iprint = argc - 1;
135  Teuchos::RCP<std::ostream> outStream;
136  Teuchos::oblackholestream bhs; // outputs nothing
137  if (iprint > 3)
138  outStream = Teuchos::rcp(&std::cout, false);
139  else
140  outStream = Teuchos::rcp(&bhs, false);
141 
142  // Save the format state of the original std::cout.
143  Teuchos::oblackholestream oldFormatState;
144  oldFormatState.copyfmt(std::cout);
145 
146  *outStream \
147  << "===============================================================================\n" \
148  << "| |\n" \
149  << "| Example: Generate Stiffness Matrix and Right Hand Side Vector for |\n" \
150  << "| Poisson Equation on Hexahedral Mesh |\n" \
151  << "| |\n" \
152  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
153  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
154  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
155  << "| |\n" \
156  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
157  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
158  << "| |\n" \
159  << "===============================================================================\n";
160 
161 
162 // ************************************ GET INPUTS **************************************
163 
164  int NX = atoi(argv[1]); // num intervals in x direction (assumed box domain, 0,1)
165  int NY = atoi(argv[2]); // num intervals in y direction (assumed box domain, 0,1)
166  int NZ = atoi(argv[3]); // num intervals in z direction (assumed box domain, 0,1)
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 <<
185  std::setw(5) << NZ << "\n\n";
186 
187  // Print mesh information
188  int numElems = NX*NY*NZ;
189  int numNodes = (NX+1)*(NY+1)*(NZ+1);
190  *outStream << " Number of Elements: " << numElems << " \n";
191  *outStream << " Number of Nodes: " << numNodes << " \n\n";
192 
193  // Cube
194  double leftX = 0.0, rightX = 1.0;
195  double leftY = 0.0, rightY = 1.0;
196  double leftZ = 0.0, rightZ = 1.0;
197 
198  // Mesh spacing
199  double hx = (rightX-leftX)/((double)NX);
200  double hy = (rightY-leftY)/((double)NY);
201  double hz = (rightZ-leftZ)/((double)NZ);
202 
203  // Get nodal coordinates
204  FieldContainer<double> nodeCoord(numNodes, spaceDim);
205  FieldContainer<int> nodeOnBoundary(numNodes);
206  int inode = 0;
207  for (int k=0; k<NZ+1; k++) {
208  for (int j=0; j<NY+1; j++) {
209  for (int i=0; i<NX+1; i++) {
210  nodeCoord(inode,0) = leftX + (double)i*hx;
211  nodeCoord(inode,1) = leftY + (double)j*hy;
212  nodeCoord(inode,2) = leftZ + (double)k*hz;
213  if (k==0 || j==0 || i==0 || k==NZ || j==NY || i==NX){
214  nodeOnBoundary(inode)=1;
215  }
216  else {
217  nodeOnBoundary(inode)=0;
218  }
219  inode++;
220  }
221  }
222  }
223 #define DUMP_DATA
224 #ifdef DUMP_DATA
225  // Print nodal coords
226  ofstream fcoordout("coords.dat");
227  for (int i=0; i<numNodes; i++) {
228  fcoordout << nodeCoord(i,0) <<" ";
229  fcoordout << nodeCoord(i,1) <<" ";
230  fcoordout << nodeCoord(i,2) <<"\n";
231  }
232  fcoordout.close();
233 #endif
234 
235 
236  // Element to Node map
237  FieldContainer<int> elemToNode(numElems, numNodesPerElem);
238  int ielem = 0;
239  for (int k=0; k<NZ; k++) {
240  for (int j=0; j<NY; j++) {
241  for (int i=0; i<NX; i++) {
242  elemToNode(ielem,0) = (NY + 1)*(NX + 1)*k + (NX + 1)*j + i;
243  elemToNode(ielem,1) = (NY + 1)*(NX + 1)*k + (NX + 1)*j + i + 1;
244  elemToNode(ielem,2) = (NY + 1)*(NX + 1)*k + (NX + 1)*(j + 1) + i + 1;
245  elemToNode(ielem,3) = (NY + 1)*(NX + 1)*k + (NX + 1)*(j + 1) + i;
246  elemToNode(ielem,4) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*j + i;
247  elemToNode(ielem,5) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*j + i + 1;
248  elemToNode(ielem,6) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*(j + 1) + i + 1;
249  elemToNode(ielem,7) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*(j + 1) + i;
250  ielem++;
251  }
252  }
253  }
254 #ifdef DUMP_DATA
255  // Output connectivity
256  ofstream fe2nout("elem2node.dat");
257  for (int k=0; k<NZ; k++) {
258  for (int j=0; j<NY; j++) {
259  for (int i=0; i<NX; i++) {
260  int ielem = i + j * NX + k * NX * NY;
261  for (int m=0; m<numNodesPerElem; m++){
262  fe2nout << elemToNode(ielem,m) <<" ";
263  }
264  fe2nout <<"\n";
265  }
266  }
267  }
268  fe2nout.close();
269 #endif
270 
271 
272 // ************************************ CUBATURE **************************************
273 
274  *outStream << "Getting cubature ... \n\n";
275 
276  // Get numerical integration points and weights
277  DefaultCubatureFactory<double> cubFactory;
278  int cubDegree = 2;
279  Teuchos::RCP<Cubature<double> > hexCub = cubFactory.create(hex_8, cubDegree);
280 
281  int cubDim = hexCub->getDimension();
282  int numCubPoints = hexCub->getNumPoints();
283 
284  FieldContainer<double> cubPoints(numCubPoints, cubDim);
285  FieldContainer<double> cubWeights(numCubPoints);
286 
287  hexCub->getCubature(cubPoints, cubWeights);
288 
289 
290 // ************************************** BASIS ***************************************
291 
292  *outStream << "Getting basis ... \n\n";
293 
294  // Define basis
296  int numFieldsG = hexHGradBasis.getCardinality();
297  FieldContainer<double> hexGVals(numFieldsG, numCubPoints);
298  FieldContainer<double> hexGrads(numFieldsG, numCubPoints, spaceDim);
299 
300  // Evaluate basis values and gradients at cubature points
301  hexHGradBasis.getValues(hexGVals, cubPoints, OPERATOR_VALUE);
302  hexHGradBasis.getValues(hexGrads, cubPoints, OPERATOR_GRAD);
303 
304 
305 // ******** LOOP OVER ELEMENTS TO CREATE LOCAL STIFFNESS MATRIX *************
306 
307  *outStream << "Building stiffness matrix and right hand side ... \n\n";
308 
309  // Settings and data structures for mass and stiffness matrices
311  typedef FunctionSpaceTools fst;
312  int numCells = 1;
313 
314  // Container for nodes
315  FieldContainer<double> hexNodes(numCells, numNodesPerElem, spaceDim);
316  // Containers for Jacobian
317  FieldContainer<double> hexJacobian(numCells, numCubPoints, spaceDim, spaceDim);
318  FieldContainer<double> hexJacobInv(numCells, numCubPoints, spaceDim, spaceDim);
319  FieldContainer<double> hexJacobDet(numCells, numCubPoints);
320  // Containers for element HGRAD stiffness matrix
321  FieldContainer<double> localStiffMatrix(numCells, numFieldsG, numFieldsG);
322  FieldContainer<double> weightedMeasure(numCells, numCubPoints);
323  FieldContainer<double> hexGradsTransformed(numCells, numFieldsG, numCubPoints, spaceDim);
324  FieldContainer<double> hexGradsTransformedWeighted(numCells, numFieldsG, numCubPoints, spaceDim);
325  // Containers for right hand side vectors
326  FieldContainer<double> rhsData(numCells, numCubPoints);
327  FieldContainer<double> localRHS(numCells, numFieldsG);
328  FieldContainer<double> hexGValsTransformed(numCells, numFieldsG, numCubPoints);
329  FieldContainer<double> hexGValsTransformedWeighted(numCells, numFieldsG, numCubPoints);
330  // Container for cubature points in physical space
331  FieldContainer<double> physCubPoints(numCells, numCubPoints, cubDim);
332 
333  // Global arrays in Epetra format
334  Epetra_SerialComm Comm;
335  Epetra_Map globalMapG(numNodes, 0, Comm);
336  Epetra_FECrsMatrix StiffMatrix(Copy, globalMapG, numFieldsG);
337  Epetra_FEVector rhs(globalMapG);
338 
339  // *** Element loop ***
340  for (int k=0; k<numElems; k++) {
341 
342  // Physical cell coordinates
343  for (int i=0; i<numNodesPerElem; i++) {
344  hexNodes(0,i,0) = nodeCoord(elemToNode(k,i),0);
345  hexNodes(0,i,1) = nodeCoord(elemToNode(k,i),1);
346  hexNodes(0,i,2) = nodeCoord(elemToNode(k,i),2);
347  }
348 
349  // Compute cell Jacobians, their inverses and their determinants
350  CellTools::setJacobian(hexJacobian, cubPoints, hexNodes, hex_8);
351  CellTools::setJacobianInv(hexJacobInv, hexJacobian );
352  CellTools::setJacobianDet(hexJacobDet, hexJacobian );
353 
354 // ************************** Compute element HGrad stiffness matrices *******************************
355 
356  // transform to physical coordinates
357  fst::HGRADtransformGRAD<double>(hexGradsTransformed, hexJacobInv, hexGrads);
358 
359  // compute weighted measure
360  fst::computeCellMeasure<double>(weightedMeasure, hexJacobDet, cubWeights);
361 
362  // multiply values with weighted measure
363  fst::multiplyMeasure<double>(hexGradsTransformedWeighted,
364  weightedMeasure, hexGradsTransformed);
365 
366  // integrate to compute element stiffness matrix
367  fst::integrate<double>(localStiffMatrix,
368  hexGradsTransformed, hexGradsTransformedWeighted, COMP_BLAS);
369 
370  // assemble into global matrix
371  for (int row = 0; row < numFieldsG; row++){
372  for (int col = 0; col < numFieldsG; col++){
373  int rowIndex = elemToNode(k,row);
374  int colIndex = elemToNode(k,col);
375  double val = localStiffMatrix(0,row,col);
376  StiffMatrix.InsertGlobalValues(1, &rowIndex, 1, &colIndex, &val);
377  }
378  }
379 
380 // ******************************* Build right hand side ************************************
381 
382  // transform integration points to physical points
383  CellTools::mapToPhysicalFrame(physCubPoints, cubPoints, hexNodes, hex_8);
384 
385  // evaluate right hand side function at physical points
386  for (int nPt = 0; nPt < numCubPoints; nPt++){
387 
388  double x = physCubPoints(0,nPt,0);
389  double y = physCubPoints(0,nPt,1);
390  double z = physCubPoints(0,nPt,2);
391 
392  rhsData(0,nPt) = evalDivGradu(x, y, z);
393  }
394 
395  // transform basis values to physical coordinates
396  fst::HGRADtransformVALUE<double>(hexGValsTransformed, hexGVals);
397 
398  // multiply values with weighted measure
399  fst::multiplyMeasure<double>(hexGValsTransformedWeighted,
400  weightedMeasure, hexGValsTransformed);
401 
402  // integrate rhs term
403  fst::integrate<double>(localRHS, rhsData, hexGValsTransformedWeighted,
404  COMP_BLAS);
405 
406  // assemble into global vector
407  for (int row = 0; row < numFieldsG; row++){
408  int rowIndex = elemToNode(k,row);
409  double val = -localRHS(0,row);
410  rhs.SumIntoGlobalValues(1, &rowIndex, &val);
411  }
412 
413  } // *** end element loop ***
414 
415 
416  // Assemble global matrices
417  StiffMatrix.GlobalAssemble(); StiffMatrix.FillComplete();
418  rhs.GlobalAssemble();
419 
420 
421  // Adjust stiffness matrix and rhs based on boundary conditions
422  for (int row = 0; row<numNodes; row++){
423  if (nodeOnBoundary(row)) {
424  int rowindex = row;
425  for (int col=0; col<numNodes; col++){
426  double val = 0.0;
427  int colindex = col;
428  StiffMatrix.ReplaceGlobalValues(1, &rowindex, 1, &colindex, &val);
429  }
430  double val = 1.0;
431  StiffMatrix.ReplaceGlobalValues(1, &rowindex, 1, &rowindex, &val);
432  val = 0.0;
433  rhs.ReplaceGlobalValues(1, &rowindex, &val);
434  }
435  }
436 
437 #ifdef DUMP_DATA
438  // Dump matrices to disk
439  EpetraExt::RowMatrixToMatlabFile("stiff_matrix.dat",StiffMatrix);
440  EpetraExt::MultiVectorToMatrixMarketFile("rhs_vector.dat",rhs,0,0,false);
441 #endif
442 
443  std::cout << "End Result: TEST PASSED\n";
444 
445  // reset format state of std::cout
446  std::cout.copyfmt(oldFormatState);
447 
448  return 0;
449 }
450 
451 
452 // Calculates value of exact solution u
453  double evalu(double & x, double & y, double & z)
454  {
455  /*
456  // function1
457  double exactu = sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
458  */
459 
460  // function2
461  double exactu = sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z);
462 
463  return exactu;
464  }
465 
466 // Calculates gradient of exact solution u
467  int evalGradu(double & x, double & y, double & z, double & gradu1, double & gradu2, double & gradu3)
468  {
469  /*
470  // function 1
471  gradu1 = M_PI*cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
472  gradu2 = M_PI*sin(M_PI*x)*cos(M_PI*y)*sin(M_PI*z);
473  gradu3 = M_PI*sin(M_PI*x)*sin(M_PI*y)*cos(M_PI*z);
474  */
475 
476  // function2
477  gradu1 = (M_PI*cos(M_PI*x)+sin(M_PI*x))
478  *sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z);
479  gradu2 = (M_PI*cos(M_PI*y)+sin(M_PI*y))
480  *sin(M_PI*x)*sin(M_PI*z)*exp(x+y+z);
481  gradu3 = (M_PI*cos(M_PI*z)+sin(M_PI*z))
482  *sin(M_PI*x)*sin(M_PI*y)*exp(x+y+z);
483 
484  return 0;
485  }
486 
487 // Calculates Laplacian of exact solution u
488  double evalDivGradu(double & x, double & y, double & z)
489  {
490  /*
491  // function 1
492  double divGradu = -3.0*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
493  */
494 
495  // function 2
496  double divGradu = -3.0*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z)
497  + 2.0*M_PI*cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z)
498  + 2.0*M_PI*cos(M_PI*y)*sin(M_PI*x)*sin(M_PI*z)*exp(x+y+z)
499  + 2.0*M_PI*cos(M_PI*z)*sin(M_PI*x)*sin(M_PI*y)*exp(x+y+z)
500  + 3.0*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z);
501 
502  return divGradu;
503  }
504 
virtual int getCardinality() const
Returns cardinality of the basis.
Header file for the Intrepid::CellTools class.
Header file for utility class to provide multidimensional containers.
Header file for utility class to provide array tools, such as tensor contractions, etc.
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.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Hexahedron cell...
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.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
A stateless class for operations on cell data. Provides methods for: