Intrepid
example_17.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
81// Intrepid includes
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_Time.hpp"
102#include "Teuchos_SerialDenseMatrix.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;
112using Teuchos::rcp;
113
114int main(int argc, char *argv[]) {
115
116 //Check number of arguments
117 if (argc < 4) {
118 std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
119 std::cout <<"Usage:\n\n";
120 std::cout <<" ./Intrepid_example_Drivers_Example_14.exe deg NX NY NZ verbose\n\n";
121 std::cout <<" where \n";
122 std::cout <<" int deg - polynomial degree to be used (assumed >= 1) \n";
123 std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n";
124 std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n";
125 std::cout <<" int NZ - num intervals in y direction (assumed box domain, 0,1) \n";
126 std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
127 exit(1);
128 }
129
130 // This little trick lets us print to std::cout only if
131 // a (dummy) command-line argument is provided.
132 int iprint = argc - 1;
133 Teuchos::RCP<std::ostream> outStream;
134 Teuchos::oblackholestream bhs; // outputs nothing
135 if (iprint > 2)
136 outStream = Teuchos::rcp(&std::cout, false);
137 else
138 outStream = Teuchos::rcp(&bhs, false);
139
140 // Save the format state of the original std::cout.
141 Teuchos::oblackholestream oldFormatState;
142 oldFormatState.copyfmt(std::cout);
143
144 *outStream \
145 << "===============================================================================\n" \
146 << "| |\n" \
147 << "| Example: Apply Stiffness Matrix for |\n" \
148 << "| Poisson Operator on Hexahedral Mesh with BLAS |\n" \
149 << "| |\n" \
150 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
151 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
152 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
153 << "| |\n" \
154 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
155 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
156 << "| |\n" \
157 << "===============================================================================\n";
158
159
160 // ************************************ GET INPUTS **************************************
161
162 int deg = atoi(argv[1]); // polynomial degree to use
163 int NX = atoi(argv[2]); // num intervals in x direction (assumed box domain, 0,1)
164 int NY = atoi(argv[3]); // num intervals in y direction (assumed box domain, 0,1)
165 int NZ = atoi(argv[4]); // num intervals in y direction (assumed box domain, 0,1)
166
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 << std::setw(5) << NZ << "\n\n";
185
186 // Print mesh information
187 int numElems = NX*NY*NZ;
188 int numNodes = (NX+1)*(NY+1)*(NZ+1);
189 *outStream << " Number of Elements: " << numElems << " \n";
190 *outStream << " Number of Nodes: " << numNodes << " \n\n";
191
192 // Cube
193 double leftX = 0.0, rightX = 1.0;
194 double leftY = 0.0, rightY = 1.0;
195 double leftZ = 0.0, rightZ = 1.0;
196
197 // Mesh spacing
198 double hx = (rightX-leftX)/((double)NX);
199 double hy = (rightY-leftY)/((double)NY);
200 double hz = (rightZ-leftZ)/((double)NZ);
201
202 // Get nodal coordinates
203 FieldContainer<double> nodeCoord(numNodes, spaceDim);
204 FieldContainer<int> nodeOnBoundary(numNodes);
205 int inode = 0;
206 for (int k=0; k<NZ+1; k++)
207 {
208 for (int j=0; j<NY+1; j++)
209 {
210 for (int i=0; i<NX+1; i++)
211 {
212 nodeCoord(inode,0) = leftX + (double)i*hx;
213 nodeCoord(inode,1) = leftY + (double)j*hy;
214 nodeCoord(inode,2) = leftZ + (double)k*hz;
215 if (k==0 || k==NZ || j==0 || i==0 || j==NY || i==NX)
216 {
217 nodeOnBoundary(inode)=1;
218 }
219 else
220 {
221 nodeOnBoundary(inode)=0;
222 }
223 inode++;
224 }
225 }
226 }
227//#define DUMP_DATA
228#ifdef DUMP_DATA
229 // Print nodal coords
230 ofstream fcoordout("coords.dat");
231 for (int i=0; i<numNodes; i++) {
232 fcoordout << nodeCoord(i,0) <<" ";
233 fcoordout << nodeCoord(i,1) <<" ";
234 fcoordout << nodeCoord(i,2) <<"\n";
235 }
236 fcoordout.close();
237#endif
238
239
240
241 // ********************************* CUBATURE AND BASIS ***********************************
242 *outStream << "Getting cubature and basis ... \n\n";
243
244 // Get numerical integration points and weights
245 // I only need this on the line since I'm doing tensor products
246 Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > glcub
247 = Teuchos::rcp(new CubaturePolylib<double,FieldContainer<double>,FieldContainer<double> >(2*deg-1,PL_GAUSS_LOBATTO) );
248
249 const int numCubPoints1D = glcub->getNumPoints();
250
251 FieldContainer<double> cubPoints1D(numCubPoints1D, 1);
252 FieldContainer<double> cubWeights1D(numCubPoints1D);
253
254 glcub->getCubature(cubPoints1D,cubWeights1D);
255
256 std::vector<Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > >
257 cub_to_tensor;
258 cub_to_tensor.push_back( glcub );
259 cub_to_tensor.push_back( glcub );
260 cub_to_tensor.push_back( glcub );
261
262 CubatureTensor<double,FieldContainer<double>,FieldContainer<double> > cubhex( cub_to_tensor );
263
264 Basis_HGRAD_HEX_Cn_FEM<double, FieldContainer<double> > hexBasis( deg , POINTTYPE_SPECTRAL );
265 const int numFields = hexBasis.getCardinality();
266
267 // get the bases tabulated at the quadrature points, dimension-by-dimension
268 const int numCubPoints = cubhex.getNumPoints();
269 FieldContainer<double> cubPoints3D( numCubPoints , spaceDim );
270 FieldContainer<double> cubWeights3D( numCubPoints );
271 FieldContainer<double> basisGrads( numFields , numCubPoints , spaceDim );
272 cubhex.getCubature( cubPoints3D , cubWeights3D );
273 hexBasis.getValues( basisGrads , cubPoints3D , OPERATOR_GRAD );
274
275
276 FieldContainer<int> elemToNode(numElems, numNodesPerElem);
277 int ielem = 0;
278 for (int k=0; k<NZ; k++)
279 {
280 for (int j=0; j<NY; j++)
281 {
282 for (int i=0; i<NX; i++)
283 {
284 elemToNode(ielem,0) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
285 elemToNode(ielem,1) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
286 elemToNode(ielem,2) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
287 elemToNode(ielem,3) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
288 elemToNode(ielem,4) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
289 elemToNode(ielem,5) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
290 elemToNode(ielem,6) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
291 elemToNode(ielem,7) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
292 ielem++;
293 }
294 }
295 }
296#ifdef DUMP_DATA
297 // Output connectivity
298 ofstream fe2nout("elem2node.dat");
299 for (int k=0;k<NZ;k++)
300 {
301 for (int j=0; j<NY; j++)
302 {
303 for (int i=0; i<NX; i++)
304 {
305 ielem = i + j * NX + k * NY * NY;
306 for (int m=0; m<numNodesPerElem; m++)
307 {
308 fe2nout << elemToNode(ielem,m) <<" ";
309 }
310 fe2nout <<"\n";
311 }
312 }
313 }
314 fe2nout.close();
315#endif
316
317
318 // ********************************* 3-D LOCAL-TO-GLOBAL MAPPING *******************************
319 FieldContainer<int> ltgMapping(numElems,numFields);
320 const int numDOF = (NX*deg+1)*(NY*deg+1)*(NZ*deg+1);
321 ielem=0;
322 for (int k=0;k<NZ;k++)
323 {
324 for (int j=0;j<NY;j++)
325 {
326 for (int i=0;i<NX;i++)
327 {
328 const int start = k * ( NY * deg + 1 ) * ( NX * deg + 1 ) + j * ( NX * deg + 1 ) + i * deg;
329 // loop over local dof on this cell
330 int local_dof_cur=0;
331 for (int kloc=0;kloc<=deg;kloc++)
332 {
333 for (int jloc=0;jloc<=deg;jloc++)
334 {
335 for (int iloc=0;iloc<=deg;iloc++)
336 {
337 ltgMapping(ielem,local_dof_cur) = start
338 + kloc * ( NX * deg + 1 ) * ( NY * deg + 1 )
339 + jloc * ( NX * deg + 1 )
340 + iloc;
341 local_dof_cur++;
342 }
343 }
344 }
345 ielem++;
346 }
347 }
348 }
349#ifdef DUMP_DATA
350 // Output ltg mapping
351 ielem = 0;
352 ofstream ltgout("ltg.dat");
353 for (int k=0;k<NZ;k++)
354 {
355 for (int j=0; j<NY; j++)
356 {
357 for (int i=0; i<NX; i++)
358 {
359 ielem = i + j * NX + k * NX * NY;
360 for (int m=0; m<numFields; m++)
361 {
362 ltgout << ltgMapping(ielem,m) <<" ";
363 }
364 ltgout <<"\n";
365 }
366 }
367 }
368 ltgout.close();
369#endif
370
371 // ********** DECLARE GLOBAL OBJECTS *************
372 Epetra_SerialComm Comm;
373 Epetra_Map globalMapG(numDOF, 0, Comm);
374 Epetra_FEVector u(globalMapG); u.Random();
375 Epetra_FEVector Ku(globalMapG);
376
377 // ************* For Jacobians **********************
378 FieldContainer<double> cellVertices(numElems,numNodesPerElem,spaceDim);
379 FieldContainer<double> cellJacobian(numElems,numCubPoints,spaceDim,spaceDim);
380 FieldContainer<double> cellJacobInv(numElems,numCubPoints,spaceDim,spaceDim);
381 FieldContainer<double> cellJacobDet(numElems,numCubPoints);
382
383
384 // get vertices of cells (for computing Jacobians)
385 for (int i=0;i<numElems;i++)
386 {
387 for (int j=0;j<numNodesPerElem;j++)
388 {
389 const int nodeCur = elemToNode(i,j);
390 for (int k=0;k<spaceDim;k++)
391 {
392 cellVertices(i,j,k) = nodeCoord(nodeCur,k);
393 }
394 }
395 }
396
397 // jacobian evaluation
398 CellTools<double>::setJacobian(cellJacobian,cubPoints3D,cellVertices,hex_8);
399 CellTools<double>::setJacobianInv(cellJacobInv, cellJacobian );
400 CellTools<double>::setJacobianDet(cellJacobDet, cellJacobian );
401
402
403 // ************* MATRIX-FREE APPLICATION
404 FieldContainer<double> uScattered(numElems,1,numFields);
405 FieldContainer<double> KuScattered(numElems,1,numFields);
406 FieldContainer<double> gradU(numElems,1,numFields,3);
407
408 // get views in SDM format for BLAS purposes
409 Teuchos::SerialDenseMatrix<int,double> refGradsAsMat( Teuchos::View ,
410 &basisGrads[0] ,
411 numCubPoints * spaceDim ,
412 numCubPoints * spaceDim ,
413 numFields );
414 Teuchos::SerialDenseMatrix<int,double> uScatAsMat( Teuchos::View ,
415 &uScattered[0] ,
416 numFields ,
417 numFields ,
418 numElems );
419 Teuchos::SerialDenseMatrix<int,double> GraduScatAsMat( Teuchos::View ,
420 &gradU[0] ,
421 numCubPoints * spaceDim ,
422 numCubPoints * spaceDim ,
423 numElems );
424 Teuchos::SerialDenseMatrix<int,double> KuScatAsMat( Teuchos::View ,
425 &KuScattered[0] ,
426 numFields ,
427 numFields ,
428 numElems );
429
430
431 u.GlobalAssemble();
432
433
434
435 Ku.PutScalar(0.0);
436 Ku.GlobalAssemble();
437
438 double *uVals = u[0];
439 double *KuVals = Ku[0];
440
441 Teuchos::Time full_timer( "Time to apply operator matrix-free:" );
442 Teuchos::Time scatter_timer( "Time to scatter dof:" );
443 Teuchos::Time elementwise_timer( "Time to do elementwise computation:" );
444 Teuchos::Time grad_timer( "Time to compute gradients:" );
445 Teuchos::Time pointwise_timer( "Time to do pointwise transformations:" );
446 Teuchos::Time moment_timer( "Time to compute moments:" );
447 Teuchos::Time gather_timer( "Time to gather dof:" );
448
449 full_timer.start();
450
451 scatter_timer.start();
452 for (int k=0; k<numElems; k++)
453 {
454 for (int i=0;i<numFields;i++)
455 {
456 uScattered(k,0,i) = uVals[ltgMapping(k,i)];
457 }
458 }
459 scatter_timer.stop();
460
461 elementwise_timer.start();
462
463 grad_timer.start();
464 // reference grad per cell with BLAS
465 GraduScatAsMat.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS ,
466 1.0 , refGradsAsMat , uScatAsMat , 0.0 );
467
468
469 grad_timer.stop();
470 pointwise_timer.start();
471 // pointwise transformations of gradients
472 Intrepid::FunctionSpaceToolsInPlace::HGRADtransformGRAD<double>( gradU , cellJacobian );
473 Intrepid::FunctionSpaceToolsInPlace::HGRADtransformGRADDual<double>( gradU , cellJacobian );
474 Intrepid::FunctionSpaceToolsInPlace::multiplyMeasure<double>( gradU , cellJacobDet );
475 pointwise_timer.stop();
476 moment_timer.start();
477
478 KuScatAsMat.multiply( Teuchos::TRANS , Teuchos::NO_TRANS , 1.0 , refGradsAsMat , GraduScatAsMat , 0.0 );
479
480 moment_timer.stop();
481 elementwise_timer.stop();
482 gather_timer.start();
483 for (int k=0;k<numElems;k++)
484 {
485 for (int i=0;i<numFields;i++)
486 {
487 KuVals[ltgMapping(k,i)] += KuScattered(k,0,i);
488 }
489 }
490 gather_timer.stop();
491 full_timer.stop();
492
493 *outStream << full_timer.name() << " " << full_timer.totalElapsedTime() << " sec\n";
494 *outStream << "\t" << scatter_timer.name() << " " << scatter_timer.totalElapsedTime() << " sec\n";
495 *outStream << "\t" << elementwise_timer.name() << " " << elementwise_timer.totalElapsedTime() << " sec\n";
496 *outStream << "\t\t" << grad_timer.name() << " " << grad_timer.totalElapsedTime() << " sec\n";
497 *outStream << "\t\t" << pointwise_timer.name() << " " << pointwise_timer.totalElapsedTime() << " sec\n";
498 *outStream << "\t\t" << moment_timer.name() << " " << moment_timer.totalElapsedTime() << " sec\n";
499 *outStream << "\t" << gather_timer.name() << " " << gather_timer.totalElapsedTime() << " sec\n";
500
501
502 *outStream << "End Result: TEST PASSED\n";
503
504 // reset format state of std::cout
505 std::cout.copyfmt(oldFormatState);
506
507 return 0;
508}
509
Header file for the Intrepid::CellTools class.
Header file for the Intrepid::CubaturePolylib class.
Header file for the Intrepid::CubatureTensor class.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::FunctionSpaceToolsInPlace class.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for the Intrepid::HGRAD_HEX_Cn_FEM class.
Header file for the Intrepid::TensorProductSpaceTools 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.
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin,...