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
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
110using namespace std;
111using namespace Intrepid;
112
113int 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";
374 typedef CellTools<double> CellTools;
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 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.
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.
Defines expert-level interfaces for the evaluation of functions and operators in physical space (supp...