Intrepid
example_03.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
80// Intrepid includes
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
109using namespace std;
110using namespace Intrepid;
111
112// Functions to evaluate exact solution and derivatives
113double evalu(double & x, double & y, double & z);
114int evalGradu(double & x, double & y, double & z, double & gradu1, double & gradu2, double & gradu3);
115double evalDivGradu(double & x, double & y, double & z);
116
117int 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 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
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
Header file for utility class to provide array tools, such as tensor contractions,...
Header file for the Intrepid::CellTools class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Intrepid utilities.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Hexahedron cell.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
virtual int getCardinality() const
Returns cardinality of the basis.
A stateless class for operations on cell data. Provides methods for:
static void mapToPhysicalFrame(ArrayPhysPoint &physPoints, const ArrayRefPoint &refPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computes F, the reference-to-physical frame map.
static void setJacobianDet(ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F.
static void setJacobianInv(ArrayJacInv &jacobianInv, const ArrayJac &jacobian)
Computes the inverse of the Jacobian matrix DF of the reference-to-physical frame map F.
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.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
Defines expert-level interfaces for the evaluation of functions and operators in physical space (supp...