Intrepid
example_16.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
103// Shards includes
104#include "Shards_CellTopology.hpp"
105
106// EpetraExt includes
107#include "EpetraExt_MultiVectorOut.h"
108
109using namespace std;
110using namespace Intrepid;
111using Teuchos::rcp;
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_14.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: Apply 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
240 // ********************************* CUBATURE AND BASIS ***********************************
241 *outStream << "Getting cubature and basis ... \n\n";
242
243 // Get numerical integration points and weights
244 // I only need this on the line since I'm doing tensor products
245 Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > glcub
246 = Teuchos::rcp(new CubaturePolylib<double,FieldContainer<double>,FieldContainer<double> >(2*deg-1,PL_GAUSS_LOBATTO) );
247
248 const int numCubPoints1D = glcub->getNumPoints();
249
250 FieldContainer<double> cubPoints1D(numCubPoints1D, 1);
251 FieldContainer<double> cubWeights1D(numCubPoints1D);
252
253 glcub->getCubature(cubPoints1D,cubWeights1D);
254
255 std::vector<Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > >
256 cub_to_tensor;
257 cub_to_tensor.push_back( glcub );
258 cub_to_tensor.push_back( glcub );
259 cub_to_tensor.push_back( glcub );
260
261 Array<RCP<FieldContainer<double> > > wts_by_dim(3);
262 wts_by_dim[0] = rcp( &cubWeights1D , false ); wts_by_dim[1] = wts_by_dim[0]; wts_by_dim[2] = wts_by_dim[1];
263
265
266 Basis_HGRAD_HEX_Cn_FEM<double, FieldContainer<double> > hexBasis( deg , POINTTYPE_SPECTRAL );
267
268 Array< Array< RCP< Basis< double , FieldContainer<double> > > > > &bases = hexBasis.getBases();
269
270 // get the bases tabulated at the quadrature points, dimension-by-dimension
271
272 Array< RCP< FieldContainer<double> > > basisVals( 3 );
273 FieldContainer<double> bvals1D( bases[0][0]->getCardinality() , numCubPoints1D );
274 bases[0][0]->getValues( bvals1D , cubPoints1D , OPERATOR_VALUE );
275 basisVals[0] = rcp( &bvals1D , false ); basisVals[1] = basisVals[0]; basisVals[2] = basisVals[0];
276
277 Array< RCP< FieldContainer<double> > > basisDVals( 3 );
278 FieldContainer<double> bdvals1D( bases[0][0]->getCardinality() , numCubPoints1D , 1);
279 bases[0][0]->getValues( bdvals1D , cubPoints1D , OPERATOR_D1 );
280 basisDVals[0] = rcp( &bdvals1D , false ); basisDVals[1] = basisDVals[0]; basisDVals[2] = basisDVals[0];
281
282
283 const int numCubPoints = cubhex.getNumPoints();
284 FieldContainer<double> cubPoints3D( numCubPoints , 3 );
285 FieldContainer<double> cubWeights3D( numCubPoints );
286 cubhex.getCubature( cubPoints3D , cubWeights3D );
287
288
289 FieldContainer<int> elemToNode(numElems, numNodesPerElem);
290 int ielem = 0;
291 for (int k=0; k<NZ; k++)
292 {
293 for (int j=0; j<NY; j++)
294 {
295 for (int i=0; i<NX; i++)
296 {
297 elemToNode(ielem,0) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
298 elemToNode(ielem,1) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
299 elemToNode(ielem,2) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
300 elemToNode(ielem,3) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
301 elemToNode(ielem,4) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
302 elemToNode(ielem,5) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
303 elemToNode(ielem,6) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
304 elemToNode(ielem,7) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
305 ielem++;
306 }
307 }
308 }
309#ifdef DUMP_DATA
310 // Output connectivity
311 ofstream fe2nout("elem2node.dat");
312 for (int k=0;k<NZ;k++)
313 {
314 for (int j=0; j<NY; j++)
315 {
316 for (int i=0; i<NX; i++)
317 {
318 ielem = i + j * NX + k * NY * NY;
319 for (int m=0; m<numNodesPerElem; m++)
320 {
321 fe2nout << elemToNode(ielem,m) <<" ";
322 }
323 fe2nout <<"\n";
324 }
325 }
326 }
327 fe2nout.close();
328#endif
329
330
331 // ********************************* 3-D LOCAL-TO-GLOBAL MAPPING *******************************
332 FieldContainer<int> ltgMapping(numElems,hexBasis.getCardinality());
333 const int numDOF = (NX*deg+1)*(NY*deg+1)*(NZ*deg+1);
334 ielem=0;
335 for (int k=0;k<NZ;k++)
336 {
337 for (int j=0;j<NY;j++)
338 {
339 for (int i=0;i<NX;i++)
340 {
341 const int start = k * ( NY * deg + 1 ) * ( NX * deg + 1 ) + j * ( NX * deg + 1 ) + i * deg;
342 // loop over local dof on this cell
343 int local_dof_cur=0;
344 for (int kloc=0;kloc<=deg;kloc++)
345 {
346 for (int jloc=0;jloc<=deg;jloc++)
347 {
348 for (int iloc=0;iloc<=deg;iloc++)
349 {
350 ltgMapping(ielem,local_dof_cur) = start
351 + kloc * ( NX * deg + 1 ) * ( NY * deg + 1 )
352 + jloc * ( NX * deg + 1 )
353 + iloc;
354 local_dof_cur++;
355 }
356 }
357 }
358 ielem++;
359 }
360 }
361 }
362#ifdef DUMP_DATA
363 // Output ltg mapping
364 ielem = 0;
365 ofstream ltgout("ltg.dat");
366 for (int k=0;k<NZ;k++)
367 {
368 for (int j=0; j<NY; j++)
369 {
370 for (int i=0; i<NX; i++)
371 {
372 ielem = i + j * NX + k * NX * NY;
373 for (int m=0; m<hexBasis.getCardinality(); m++)
374 {
375 ltgout << ltgMapping(ielem,m) <<" ";
376 }
377 ltgout <<"\n";
378 }
379 }
380 }
381 ltgout.close();
382#endif
383
384 // ********** DECLARE GLOBAL OBJECTS *************
385 Epetra_SerialComm Comm;
386 Epetra_Map globalMapG(numDOF, 0, Comm);
387 Epetra_FEVector u(globalMapG); u.Random();
388 Epetra_FEVector Ku(globalMapG);
389
390 // ************* For Jacobians **********************
391 FieldContainer<double> cellVertices(numElems,numNodesPerElem,spaceDim);
392 FieldContainer<double> cellJacobian(numElems,numCubPoints,spaceDim,spaceDim);
393 FieldContainer<double> cellJacobInv(numElems,numCubPoints,spaceDim,spaceDim);
394 FieldContainer<double> cellJacobDet(numElems,numCubPoints);
395
396
397 // get vertices of cells (for computing Jacobians)
398 for (int i=0;i<numElems;i++)
399 {
400 for (int j=0;j<numNodesPerElem;j++)
401 {
402 const int nodeCur = elemToNode(i,j);
403 for (int k=0;k<spaceDim;k++)
404 {
405 cellVertices(i,j,k) = nodeCoord(nodeCur,k);
406 }
407 }
408 }
409
410 // jacobian evaluation
411 CellTools<double>::setJacobian(cellJacobian,cubPoints3D,cellVertices,hex_8);
412 CellTools<double>::setJacobianInv(cellJacobInv, cellJacobian );
413 CellTools<double>::setJacobianDet(cellJacobDet, cellJacobian );
414
415
416 // ************* MATRIX-FREE APPLICATION
417 FieldContainer<double> uScattered(numElems,1,hexBasis.getCardinality());
418 FieldContainer<double> KuScattered(numElems,1,hexBasis.getCardinality());
419 FieldContainer<double> gradU(numElems,1,hexBasis.getCardinality(),3);
420
421 u.GlobalAssemble();
422
423
424
425 Ku.PutScalar(0.0);
426 Ku.GlobalAssemble();
427
428 double *uVals = u[0];
429 double *KuVals = Ku[0];
430
431 Teuchos::Time full_timer( "Time to apply operator matrix-free:" );
432 Teuchos::Time scatter_timer( "Time to scatter dof:" );
433 Teuchos::Time elementwise_timer( "Time to do elementwise computation:" );
434 Teuchos::Time grad_timer( "Time to compute gradients:" );
435 Teuchos::Time pointwise_timer( "Time to do pointwise transformations:" );
436 Teuchos::Time moment_timer( "Time to compute moments:" );
437 Teuchos::Time gather_timer( "Time to gather dof:" );
438
439 full_timer.start();
440
441 scatter_timer.start();
442 for (int k=0; k<numElems; k++)
443 {
444 for (int i=0;i<hexBasis.getCardinality();i++)
445 {
446 uScattered(k,0,i) = uVals[ltgMapping(k,i)];
447 }
448 }
449 scatter_timer.stop();
450
451 elementwise_timer.start();
452
453 grad_timer.start();
454 Intrepid::TensorProductSpaceTools::evaluateGradient<double>( gradU , uScattered ,basisVals , basisDVals );
455 grad_timer.stop();
456 pointwise_timer.start();
457 Intrepid::FunctionSpaceToolsInPlace::HGRADtransformGRAD<double>( gradU , cellJacobian );
458 Intrepid::FunctionSpaceToolsInPlace::HGRADtransformGRADDual<double>( gradU , cellJacobian );
459 Intrepid::FunctionSpaceToolsInPlace::multiplyMeasure<double>( gradU , cellJacobDet );
460 pointwise_timer.stop();
461 moment_timer.start();
462 Intrepid::TensorProductSpaceTools::momentsGrad<double>( KuScattered , gradU , basisVals , basisDVals , wts_by_dim );
463 moment_timer.stop();
464 elementwise_timer.stop();
465 gather_timer.start();
466 for (int k=0;k<numElems;k++)
467 {
468 for (int i=0;i<hexBasis.getCardinality();i++)
469 {
470 KuVals[ltgMapping(k,i)] += KuScattered(k,0,i);
471 }
472 }
473 gather_timer.stop();
474 full_timer.stop();
475
476 *outStream << full_timer.name() << " " << full_timer.totalElapsedTime() << " sec\n";
477 *outStream << "\t" << scatter_timer.name() << " " << scatter_timer.totalElapsedTime() << " sec\n";
478 *outStream << "\t" << elementwise_timer.name() << " " << elementwise_timer.totalElapsedTime() << " sec\n";
479 *outStream << "\t\t" << grad_timer.name() << " " << grad_timer.totalElapsedTime() << " sec\n";
480 *outStream << "\t\t" << pointwise_timer.name() << " " << pointwise_timer.totalElapsedTime() << " sec\n";
481 *outStream << "\t\t" << moment_timer.name() << " " << moment_timer.totalElapsedTime() << " sec\n";
482 *outStream << "\t" << gather_timer.name() << " " << gather_timer.totalElapsedTime() << " sec\n";
483
484
485 *outStream << "End Result: TEST PASSED\n";
486
487 // reset format state of std::cout
488 std::cout.copyfmt(oldFormatState);
489
490 return 0;
491}
492
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.
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.
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin,...
Defines tensor-product cubature (integration) rules in Intrepid.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....