Intrepid
example_14.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
87//#include "Intrepid_ArrayTools.hpp"
88#include "Intrepid_HGRAD_LINE_Cn_FEM.hpp"
89//#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_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 // 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 // ********************************* 1-D CUBATURE AND BASIS ***********************************
283 *outStream << "Getting cubature and basis ... \n\n";
284
285 // Get numerical integration points and weights
286 // I only need this on the line since I'm doing tensor products
287 Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > glcub
288 = Teuchos::rcp(new CubaturePolylib<double,FieldContainer<double>,FieldContainer<double> >(2*deg-1,PL_GAUSS_LOBATTO) );
289
290 const int numCubPoints = glcub->getNumPoints();
291
292 FieldContainer<double> cubPoints1D(numCubPoints, 1);
293 FieldContainer<double> cubWeights1D(numCubPoints);
294
295 glcub->getCubature(cubPoints1D,cubWeights1D);
296 // Define basis: I only need this on the line also
297 Basis_HGRAD_LINE_Cn_FEM<double, FieldContainer<double> > lineHGradBasis(deg,POINTTYPE_SPECTRAL);
298 int numLineFieldsG = lineHGradBasis.getCardinality();
299 FieldContainer<double> lineGrads(numLineFieldsG, numCubPoints, 1);
300
301 // Evaluate basis values and gradients at cubature points
302 lineHGradBasis.getValues(lineGrads, cubPoints1D, OPERATOR_GRAD);
303
304
305 // ********************************* 3-D LOCAL-TO-GLOBAL MAPPING *******************************
306 FieldContainer<int> ltgMapping(numElems,numLineFieldsG*numLineFieldsG*numLineFieldsG);
307 const int numDOF = (NX*deg+1)*(NY*deg+1)*(NZ*deg+1);
308 ielem=0;
309 for (int k=0;k<NZ;k++)
310 {
311 for (int j=0;j<NY;j++)
312 {
313 for (int i=0;i<NX;i++)
314 {
315 const int start = k * ( NY * deg + 1 ) * ( NX * deg + 1 ) + j * ( NX * deg + 1 ) + i * deg;
316 // loop over local dof on this cell
317 int local_dof_cur=0;
318 for (int kloc=0;kloc<=deg;kloc++)
319 {
320 for (int jloc=0;jloc<=deg;jloc++)
321 {
322 for (int iloc=0;iloc<=deg;iloc++)
323 {
324 ltgMapping(ielem,local_dof_cur) = start
325 + kloc * ( NX * deg + 1 ) * ( NY * deg + 1 )
326 + jloc * ( NX * deg + 1 )
327 + iloc;
328 local_dof_cur++;
329 }
330 }
331 }
332 ielem++;
333 }
334 }
335 }
336#ifdef DUMP_DATA
337 // Output ltg mapping
338 ofstream ltgout("ltg.dat");
339 for (int k=0;k<NZ;k++)
340 {
341 for (int j=0; j<NY; j++)
342 {
343 for (int i=0; i<NX; i++)
344 {
345 ielem = i + j * NX + k * NX * NY;
346 for (int m=0; m<numLineFieldsG*numLineFieldsG*numLineFieldsG; m++)
347 {
348 ltgout << ltgMapping(ielem,m) <<" ";
349 }
350 ltgout <<"\n";
351 }
352 }
353 }
354 ltgout.close();
355#endif
356
357 // ********** DECLARE GLOBAL OBJECTS *************
358 Epetra_SerialComm Comm;
359 Epetra_Map globalMapG(numDOF, 0, Comm);
360 Epetra_FEVector u(globalMapG); u.Random();
361 Epetra_FEVector Ku(globalMapG);
362
363
364 // ************* MATRIX-FREE APPLICATION
365 FieldContainer<double> uScattered(numElems,numLineFieldsG*numLineFieldsG*numLineFieldsG);
366 FieldContainer<double> KuScattered(numElems,numLineFieldsG*numLineFieldsG*numLineFieldsG);
367
368 u.GlobalAssemble();
369
370 Epetra_Time multTimer(Comm);
371 Teuchos::BLAS<int,double> blas;
372 Ku.PutScalar(0.0);
373 Ku.GlobalAssemble();
374
375 double *uVals = u[0];
376 double *KuVals = Ku[0];
377
378 Epetra_Time scatterTimer(Comm);
379 std::cout << "Scattering\n";
380 // Scatter
381 for (int k=0; k<numElems; k++)
382 {
383 for (int i=0;i<numLineFieldsG*numLineFieldsG*numLineFieldsG;i++)
384 {
385 uScattered(k,i) = uVals[ltgMapping(k,i)];
386 }
387 }
388
389
390 const double scatterTime = scatterTimer.ElapsedTime();
391
392
393
394 FieldContainer<double> Du(numLineFieldsG,numLineFieldsG,numLineFieldsG);
395
396 Epetra_Time localAppTimer(Comm);
397
398 uScattered.resize(numElems,numLineFieldsG,numLineFieldsG,numLineFieldsG);
399
400
401 int cur;
402 double hcur;
403
404 for (ielem=0;ielem<numElems;ielem++)
405 {
406 // X-COMPONENT OF ELEMENT LAPLACIAN
407
408 // zero out derivative
409 for (int k=0;k<numLineFieldsG;k++)
410 {
411 for (int j=0;j<numLineFieldsG;j++)
412 {
413 for (int i=0;i<numLineFieldsG;i++)
414 {
415 Du(k,j,i) = 0.0;
416 }
417 }
418 }
419
420
421 // compute x derivative
422 // this could be a very simple dgemm call
423 for (int k=0;k<numLineFieldsG;k++)
424 {
425 for (int j=0;j<numLineFieldsG;j++)
426 {
427 for (int i=0;i<numLineFieldsG;i++)
428 {
429 for (int q=0;q<numLineFieldsG;q++)
430 {
431 Du(k,j,i) += uScattered(ielem,k,j,i) * lineGrads(i,q,0);
432 }
433 }
434 }
435 }
436
437 // integration loop for x derivative term
438 cur = 0;
439 hcur = hy * hz / hx;
440 for (int k=0;k<numLineFieldsG;k++)
441 {
442 const double wt1 = hcur * cubWeights1D(k);
443 for (int j=0;j<numLineFieldsG;j++)
444 {
445 const double wtstuff = wt1 * cubWeights1D(j);
446 for (int i=0;i<numLineFieldsG;i++)
447 {
448 for (int q=0;q<numLineFieldsG;q++)
449 {
450 KuScattered(ielem,cur) += wtstuff
451 * cubWeights1D(q) * Du(k,j,q) * lineGrads(i,q,0);
452 }
453 cur++;
454 }
455 }
456 }
457
458
459 // Y-COMPONENT OF ELEMENT LAPLACIAN
460
461 // zero out derivative
462 for (int k=0;k<numLineFieldsG;k++)
463 {
464 for (int j=0;j<numLineFieldsG;j++)
465 {
466 for (int i=0;i<numLineFieldsG;i++)
467 {
468 Du(k,j,i) = 0.0;
469 }
470 }
471 }
472
473 // compute y derivative
474 for (int k=0;k<numLineFieldsG;k++)
475 {
476 for (int j=0;j<numLineFieldsG;j++)
477 {
478 for (int i=0;i<numLineFieldsG;i++)
479 {
480 for (int q=0;q<numLineFieldsG;q++)
481 {
482 Du(k,j,i) += uScattered(ielem,k,j,i) * lineGrads(j,q,0);
483 }
484 }
485 }
486 }
487
488 // integration loop for y derivative term
489 cur = 0;
490 hcur = hx * hz / hy;
491 for (int k=0;k<numLineFieldsG;k++)
492 {
493 const double wt1 = hcur * cubWeights1D(k);
494 for (int j=0;j<numLineFieldsG;j++)
495 {
496 for (int i=0;i<numLineFieldsG;i++)
497 {
498 const double wtstuff = cubWeights1D(i) * wt1;
499 for (int q=0;q<numLineFieldsG;q++)
500 {
501 KuScattered(ielem,cur) += wtstuff * cubWeights1D(q) * Du(k,q,i) * lineGrads(j,q,0);
502 }
503 cur++;
504 }
505 }
506 }
507
508
509 // Z-COMPONENT OF ELEMENT LAPLACIAN
510
511 // zero out derivative
512 for (int k=0;k<numLineFieldsG;k++)
513 {
514 for (int j=0;j<numLineFieldsG;j++)
515 {
516 for (int i=0;i<numLineFieldsG;i++)
517 {
518 Du(k,j,i) = 0.0;
519 }
520 }
521 }
522
523 // compute z derivative
524 for (int k=0;k<numLineFieldsG;k++)
525 {
526 for (int j=0;j<numLineFieldsG;j++)
527 {
528 for (int i=0;i<numLineFieldsG;i++)
529 {
530 for (int q=0;q<numLineFieldsG;q++)
531 {
532 Du(k,j,i) += uScattered(ielem,k,j,i) * lineGrads(k,q,0);
533 }
534 }
535 }
536 }
537
538 // integration loop for z derivative term.
539 cur = 0;
540 hcur = hx * hy / hz;
541 for (int k=0;k<numLineFieldsG;k++)
542 {
543 for (int j=0;j<numLineFieldsG;j++)
544 {
545 const double wt1 = hcur * cubWeights1D(j);
546 for (int i=0;i<numLineFieldsG;i++)
547 {
548 const double wtstuff = cubWeights1D(i) * wt1;
549 for (int q=0;q<numLineFieldsG;q++)
550 {
551 KuScattered(ielem,cur) += wtstuff * cubWeights1D(q) * Du(q,j,i) * lineGrads(k,q,0);
552 }
553 cur++;
554 }
555 }
556 }
557
558 }
559
560 const double localAppTime = localAppTimer.ElapsedTime();
561
562 Epetra_Time gatherTimer(Comm);
563 // Gather
564 for (int k=0;k<numElems;k++)
565 {
566 for (int i=0;i<numLineFieldsG*numLineFieldsG*numLineFieldsG;i++)
567 {
568 KuVals[ltgMapping(k,i)] += KuScattered(k,i);
569 }
570 }
571
572 const double gatherTime = gatherTimer.ElapsedTime();
573
574
575 *outStream << "Time to scatter " << scatterTime << "\n";
576 *outStream << "Time for local application " << localAppTime << "\n";
577 *outStream << "Time to gather " << gatherTime << "\n";
578 *outStream << "Total matrix-free time " << scatterTime + localAppTime + gatherTime << "\n";
579
580
581 *outStream << "End Result: TEST PASSED\n";
582
583 // reset format state of std::cout
584 std::cout.copyfmt(oldFormatState);
585
586 return 0;
587}
588
Header file for the Intrepid::CellTools class.
Header file for the Intrepid::CubaturePolylib class.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::FunctionSpaceTools class.
Intrepid utilities.
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin,...