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
85 #include "Intrepid_CellTools.hpp"
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 
110 using namespace std;
111 using namespace Intrepid;
112 
113 int 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 
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1...
Header file for the Intrepid::CellTools class.
Header file for utility class to provide multidimensional containers.
Intrepid utilities.
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin, Aeronautics, Imperial College London) within Intrepid.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for the Intrepid::CubaturePolylib class.