Intrepid
example_15.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
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_FECrsMatrix.h"
97#include "Epetra_SerialComm.h"
98
99// Teuchos includes
100#include "Teuchos_oblackholestream.hpp"
101#include "Teuchos_RCP.hpp"
102//#include "Teuchos_BLAS.hpp"
103//#include "Teuchos_BLAS_types.hpp"
104
105// Shards includes
106#include "Shards_CellTopology.hpp"
107
108// EpetraExt includes
109#include "EpetraExt_MultiVectorOut.h"
110
111#include <vector>
112#include <map>
113
114using namespace std;
115using namespace Intrepid;
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_15.exe deg NX NY NZ verbose\n\n";
124 std::cout <<" where \n";
125 std::cout <<" int deg - polynomial degree to be used (assumed >= 1) \n";
126 std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n";
127 std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n";
128 std::cout <<" int NZ - num intervals in y direction (assumed box domain, 0,1) \n";
129 std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
130 exit(1);
131 }
132
133 // This little trick lets us print to std::cout only if
134 // a (dummy) command-line argument is provided.
135 int iprint = argc - 1;
136 Teuchos::RCP<std::ostream> outStream;
137 Teuchos::oblackholestream bhs; // outputs nothing
138 if (iprint > 2)
139 outStream = Teuchos::rcp(&std::cout, false);
140 else
141 outStream = Teuchos::rcp(&bhs, false);
142
143 // Save the format state of the original std::cout.
144 Teuchos::oblackholestream oldFormatState;
145 oldFormatState.copyfmt(std::cout);
146
147 *outStream \
148 << "===============================================================================\n" \
149 << "| |\n" \
150 << "| Example: Build Stiffness Matrix for |\n" \
151 << "| Poisson Equation on Hexahedral Mesh |\n" \
152 << "| |\n" \
153 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
154 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
155 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
156 << "| |\n" \
157 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
158 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
159 << "| |\n" \
160 << "===============================================================================\n";
161
162
163 // ************************************ GET INPUTS **************************************
164
165 int deg = atoi(argv[1]); // polynomial degree to use
166 int NX = atoi(argv[2]); // num intervals in x direction (assumed box domain, 0,1)
167 int NY = atoi(argv[3]); // num intervals in y direction (assumed box domain, 0,1)
168 int NZ = atoi(argv[4]); // num intervals in y direction (assumed box domain, 0,1)
169
170
171 // *********************************** CELL TOPOLOGY **********************************
172
173 // Get cell topology for base hexahedron
174 typedef shards::CellTopology CellTopology;
175 CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() );
176
177 // Get dimensions
178 int numNodesPerElem = hex_8.getNodeCount();
179 int spaceDim = hex_8.getDimension();
180
181 // *********************************** GENERATE MESH ************************************
182
183 *outStream << "Generating mesh ... \n\n";
184
185 *outStream << " NX" << " NY" << " NZ\n";
186 *outStream << std::setw(5) << NX <<
187 std::setw(5) << NY << std::setw(5) << NZ << "\n\n";
188
189 // Print mesh information
190 int numElems = NX*NY*NZ;
191 int numNodes = (NX+1)*(NY+1)*(NZ+1);
192 *outStream << " Number of Elements: " << numElems << " \n";
193 *outStream << " Number of Nodes: " << numNodes << " \n\n";
194
195 // Cube
196 double leftX = 0.0, rightX = 1.0;
197 double leftY = 0.0, rightY = 1.0;
198 double leftZ = 0.0, rightZ = 1.0;
199
200 // Mesh spacing
201 double hx = (rightX-leftX)/((double)NX);
202 double hy = (rightY-leftY)/((double)NY);
203 double hz = (rightZ-leftZ)/((double)NZ);
204
205 // Get nodal coordinates
206 FieldContainer<double> nodeCoord(numNodes, spaceDim);
207 FieldContainer<int> nodeOnBoundary(numNodes);
208 int inode = 0;
209 for (int k=0; k<NZ+1; k++)
210 {
211 for (int j=0; j<NY+1; j++)
212 {
213 for (int i=0; i<NX+1; i++)
214 {
215 nodeCoord(inode,0) = leftX + (double)i*hx;
216 nodeCoord(inode,1) = leftY + (double)j*hy;
217 nodeCoord(inode,2) = leftZ + (double)k*hz;
218 if (k==0 || k==NZ || j==0 || i==0 || j==NY || i==NX)
219 {
220 nodeOnBoundary(inode)=1;
221 }
222 else
223 {
224 nodeOnBoundary(inode)=0;
225 }
226 inode++;
227 }
228 }
229 }
230#define DUMP_DATA
231#ifdef DUMP_DATA
232 // Print nodal coords
233 ofstream fcoordout("coords.dat");
234 for (int i=0; i<numNodes; i++) {
235 fcoordout << nodeCoord(i,0) <<" ";
236 fcoordout << nodeCoord(i,1) <<" ";
237 fcoordout << nodeCoord(i,2) <<"\n";
238 }
239 fcoordout.close();
240#endif
241
242
243 // Element to Node map
244 // We'll keep it around, but this is only the DOFMap if you are in the lowest order case.
245 FieldContainer<int> elemToNode(numElems, numNodesPerElem);
246 int ielem = 0;
247 for (int k=0; k<NZ; k++)
248 {
249 for (int j=0; j<NY; j++)
250 {
251 for (int i=0; i<NX; i++)
252 {
253 elemToNode(ielem,0) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
254 elemToNode(ielem,1) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
255 elemToNode(ielem,2) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
256 elemToNode(ielem,3) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
257 elemToNode(ielem,4) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
258 elemToNode(ielem,5) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
259 elemToNode(ielem,6) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
260 elemToNode(ielem,7) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
261 ielem++;
262 }
263 }
264 }
265#ifdef DUMP_DATA
266 // Output connectivity
267 ofstream fe2nout("elem2node.dat");
268 for (int k=0;k<NZ;k++)
269 {
270 for (int j=0; j<NY; j++)
271 {
272 for (int i=0; i<NX; i++)
273 {
274 ielem = i + j * NX + k * NY * NY;
275 for (int m=0; m<numNodesPerElem; m++)
276 {
277 fe2nout << elemToNode(ielem,m) <<" ";
278 }
279 fe2nout <<"\n";
280 }
281 }
282 }
283 fe2nout.close();
284#endif
285
286 // ************************************ CUBATURE **************************************
287 *outStream << "Getting cubature ... \n\n";
288
289 // Get numerical integration points and weights
291 int cubDegree = 2*deg;
292 Teuchos::RCP<Cubature<double> > quadCub = cubFactory.create(hex_8, cubDegree);
293
294 int cubDim = quadCub->getDimension();
295 int numCubPoints = quadCub->getNumPoints();
296
297 FieldContainer<double> cubPoints(numCubPoints, cubDim);
298 FieldContainer<double> cubWeights(numCubPoints);
299
300 quadCub->getCubature(cubPoints, cubWeights);
301
302
303 // ************************************** BASIS ***************************************
304
305 *outStream << "Getting basis ... \n\n";
306
307 // Define basis
308 Basis_HGRAD_HEX_Cn_FEM<double, FieldContainer<double> > quadHGradBasis(deg,POINTTYPE_SPECTRAL);
309 int numFieldsG = quadHGradBasis.getCardinality();
310 FieldContainer<double> quadGVals(numFieldsG, numCubPoints);
311 FieldContainer<double> quadGrads(numFieldsG, numCubPoints, spaceDim);
312
313 // Evaluate basis values and gradients at cubature points
314 quadHGradBasis.getValues(quadGVals, cubPoints, OPERATOR_VALUE);
315 quadHGradBasis.getValues(quadGrads, cubPoints, OPERATOR_GRAD);
316
317 // create the local-global mapping
318 FieldContainer<int> ltgMapping(numElems,numFieldsG);
319 const int numDOF = (NX*deg+1)*(NY*deg+1)*(NZ*deg+1);
320 ielem=0;
321 for (int k=0;k<NZ;k++)
322 {
323 for (int j=0;j<NY;j++)
324 {
325 for (int i=0;i<NX;i++)
326 {
327 const int start = k * ( NY * deg + 1 ) * ( NX * deg + 1 ) + j * ( NX * deg + 1 ) + i * deg;
328 // loop over local dof on this cell
329 int local_dof_cur=0;
330 for (int kloc=0;kloc<=deg;kloc++)
331 {
332 for (int jloc=0;jloc<=deg;jloc++)
333 {
334 for (int iloc=0;iloc<=deg;iloc++)
335 {
336 ltgMapping(ielem,local_dof_cur) = start
337 + kloc * ( NX * deg + 1 ) * ( NY * deg + 1 )
338 + jloc * ( NX * deg + 1 )
339 + iloc;
340 local_dof_cur++;
341 }
342 }
343 }
344 ielem++;
345 }
346 }
347 }
348#ifdef DUMP_DATA
349 // Output ltg mapping
350 ofstream ltgout("ltg.dat");
351 for (int k=0;k<NZ;k++)
352 {
353 for (int j=0; j<NY; j++)
354 {
355 for (int i=0; i<NX; i++)
356 {
357 ielem = i + j * NX + k * NX * NY;
358 for (int m=0; m<numFieldsG; m++)
359 {
360 ltgout << ltgMapping(ielem,m) <<" ";
361 }
362 ltgout <<"\n";
363 }
364 }
365 }
366 ltgout.close();
367#endif
368
369 // ********** DECLARE GLOBAL OBJECTS *************
370 Epetra_SerialComm Comm;
371 Epetra_Map globalMapG(numDOF, 0, Comm);
372 Epetra_FEVector u(globalMapG); u.Random();
373 Epetra_FEVector Ku(globalMapG);
374
375 // time the instantiation
376// Epetra_Time instantiateTimer(Comm);
377// Epetra_FECrsMatrix StiffMatrix(Copy,globalMapG,8*numFieldsG);
378// const double instantiateTime = instantiateTimer.ElapsedTime();
379
380
381 // ********** CONSTRUCT AND INSERT LOCAL STIFFNESS MATRICES ***********
382 *outStream << "Building local stiffness matrices...\n\n";
384 typedef FunctionSpaceTools fst;
385 int numCells = numElems;
386
387 // vertices
388 FieldContainer<double> cellVertices(numCells,numNodesPerElem,spaceDim);
389
390 // jacobian information
391 FieldContainer<double> cellJacobian(numCells,numCubPoints,spaceDim,spaceDim);
392 FieldContainer<double> cellJacobInv(numCells,numCubPoints,spaceDim,spaceDim);
393 FieldContainer<double> cellJacobDet(numCells,numCubPoints);
394
395 // element stiffness matrices and supporting storage space
396 FieldContainer<double> localStiffMatrices(numCells, numFieldsG, numFieldsG);
397 FieldContainer<double> transformedBasisGradients(numCells,numFieldsG,numCubPoints,spaceDim);
398 FieldContainer<double> weightedTransformedBasisGradients(numCells,numFieldsG,numCubPoints,spaceDim);
399 FieldContainer<double> weightedMeasure(numCells, numCubPoints);
400
401
402 // get vertices of cells (for computing Jacobians)
403 for (int i=0;i<numElems;i++)
404 {
405 for (int j=0;j<numNodesPerElem;j++)
406 {
407 const int nodeCur = elemToNode(i,j);
408 for (int k=0;k<spaceDim;k++)
409 {
410 cellVertices(i,j,k) = nodeCoord(nodeCur,k);
411 }
412 }
413 }
414
415 Epetra_Time localConstructTimer( Comm );
416
417 // jacobian evaluation
418 CellTools::setJacobian(cellJacobian,cubPoints,cellVertices,hex_8);
419 CellTools::setJacobianInv(cellJacobInv, cellJacobian );
420 CellTools::setJacobianDet(cellJacobDet, cellJacobian );
421
422 // transform reference element gradients to each cell
423 fst::HGRADtransformGRAD<double>(transformedBasisGradients, cellJacobInv, quadGrads);
424
425 // compute weighted measure
426 fst::computeCellMeasure<double>(weightedMeasure, cellJacobDet, cubWeights);
427
428 // multiply values with weighted measure
429 fst::multiplyMeasure<double>(weightedTransformedBasisGradients,
430 weightedMeasure, transformedBasisGradients);
431
432 // integrate to compute element stiffness matrix
433 fst::integrate<double>(localStiffMatrices,
434 transformedBasisGradients, weightedTransformedBasisGradients , COMP_BLAS);
435
436 const double localConstructTime = localConstructTimer.ElapsedTime();
437
438
439 Epetra_Time insertionTimer(Comm);
440
441 vector<map<int,double> > mat(numDOF);
442
443
444
445 // *** Element loop ***
446 for (int el=0; el<numElems; el++)
447 {
448 for (int i=0;i<numFieldsG;i++) // local rows
449 {
450 const int glob_row = ltgMapping(el,i);
451 map<int,double> & cur_row = mat[glob_row];
452
453 for (int j=0;j<numFieldsG;j++) // local columns
454 {
455 const int glob_col = ltgMapping(el,j);
456 const double cur_val = localStiffMatrices(el,i,j);
457 map<int,double>::iterator it = cur_row.find( glob_col );
458 if (it != cur_row.end()) // current column already in row
459 {
460 it->second += cur_val;
461 }
462 else
463 {
464 cur_row[glob_col] = cur_val;
465 }
466 }
467 }
468 }
469 //StiffMatrix.GlobalAssemble(); StiffMatrix.FillComplete();
470 const double insertionTime = insertionTimer.ElapsedTime( );
471
472// *outStream << "Time to instantiate global stiffness matrix: " << instantiateTime << "\n";
473 *outStream << "Time to build local matrices (including Jacobian computation): "<< localConstructTime << "\n";
474 *outStream << "Time to assemble global matrix from local matrices: " << insertionTime << "\n";
475 *outStream << "Total construction time: " << localConstructTime + insertionTime << "\n";
476
477// Epetra_Time applyTimer(Comm);
478// StiffMatrix.Apply(u,Ku);
479// const double multTime = applyTimer.ElapsedTime();
480// *outStream << "Time to multiply onto a vector: " << multTime << "\n";
481
482 *outStream << "End Result: TEST PASSED\n";
483
484 // reset format state of std::cout
485 std::cout.copyfmt(oldFormatState);
486
487 return 0;
488}
489
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 the Intrepid::HGRAD_HEX_Cn_FEM class.
Intrepid utilities.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell.
A stateless class for operations on cell data. Provides methods for:
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...