Intrepid
example_15.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"
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_FECrsMatrix.h"
97 #include "Epetra_SerialComm.h"
98 
99 // Teuchos includes
100 #include "Teuchos_oblackholestream.hpp"
101 #include "Teuchos_RCP.hpp"
102 //#include "Teuchos_BLAS.hpp"
103 //#include "Teuchos_BLAS_types.hpp"
104 
105 // Shards includes
106 #include "Shards_CellTopology.hpp"
107 
108 // EpetraExt includes
109 #include "EpetraExt_MultiVectorOut.h"
110 
111 #include <vector>
112 #include <map>
113 
114 using namespace std;
115 using namespace Intrepid;
116 
117 int main(int argc, char *argv[]) {
118 
119  //Check number of arguments
120  if (argc < 4) {
121  std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
122  std::cout <<"Usage:\n\n";
123  std::cout <<" ./Intrepid_example_Drivers_Example_15.exe deg NX NY NZ verbose\n\n";
124  std::cout <<" where \n";
125  std::cout <<" int deg - polynomial degree to be used (assumed >= 1) \n";
126  std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n";
127  std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n";
128  std::cout <<" int NZ - num intervals in y direction (assumed box domain, 0,1) \n";
129  std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
130  exit(1);
131  }
132 
133  // This little trick lets us print to std::cout only if
134  // a (dummy) command-line argument is provided.
135  int iprint = argc - 1;
136  Teuchos::RCP<std::ostream> outStream;
137  Teuchos::oblackholestream bhs; // outputs nothing
138  if (iprint > 2)
139  outStream = Teuchos::rcp(&std::cout, false);
140  else
141  outStream = Teuchos::rcp(&bhs, false);
142 
143  // Save the format state of the original std::cout.
144  Teuchos::oblackholestream oldFormatState;
145  oldFormatState.copyfmt(std::cout);
146 
147  *outStream \
148  << "===============================================================================\n" \
149  << "| |\n" \
150  << "| Example: Build Stiffness Matrix for |\n" \
151  << "| Poisson Equation on Hexahedral Mesh |\n" \
152  << "| |\n" \
153  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
154  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
155  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
156  << "| |\n" \
157  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
158  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
159  << "| |\n" \
160  << "===============================================================================\n";
161 
162 
163  // ************************************ GET INPUTS **************************************
164 
165  int deg = atoi(argv[1]); // polynomial degree to use
166  int NX = atoi(argv[2]); // num intervals in x direction (assumed box domain, 0,1)
167  int NY = atoi(argv[3]); // num intervals in y direction (assumed box domain, 0,1)
168  int NZ = atoi(argv[4]); // num intervals in y direction (assumed box domain, 0,1)
169 
170 
171  // *********************************** CELL TOPOLOGY **********************************
172 
173  // Get cell topology for base hexahedron
174  typedef shards::CellTopology CellTopology;
175  CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() );
176 
177  // Get dimensions
178  int numNodesPerElem = hex_8.getNodeCount();
179  int spaceDim = hex_8.getDimension();
180 
181  // *********************************** GENERATE MESH ************************************
182 
183  *outStream << "Generating mesh ... \n\n";
184 
185  *outStream << " NX" << " NY" << " NZ\n";
186  *outStream << std::setw(5) << NX <<
187  std::setw(5) << NY << std::setw(5) << NZ << "\n\n";
188 
189  // Print mesh information
190  int numElems = NX*NY*NZ;
191  int numNodes = (NX+1)*(NY+1)*(NZ+1);
192  *outStream << " Number of Elements: " << numElems << " \n";
193  *outStream << " Number of Nodes: " << numNodes << " \n\n";
194 
195  // Cube
196  double leftX = 0.0, rightX = 1.0;
197  double leftY = 0.0, rightY = 1.0;
198  double leftZ = 0.0, rightZ = 1.0;
199 
200  // Mesh spacing
201  double hx = (rightX-leftX)/((double)NX);
202  double hy = (rightY-leftY)/((double)NY);
203  double hz = (rightZ-leftZ)/((double)NZ);
204 
205  // Get nodal coordinates
206  FieldContainer<double> nodeCoord(numNodes, spaceDim);
207  FieldContainer<int> nodeOnBoundary(numNodes);
208  int inode = 0;
209  for (int k=0; k<NZ+1; k++)
210  {
211  for (int j=0; j<NY+1; j++)
212  {
213  for (int i=0; i<NX+1; i++)
214  {
215  nodeCoord(inode,0) = leftX + (double)i*hx;
216  nodeCoord(inode,1) = leftY + (double)j*hy;
217  nodeCoord(inode,2) = leftZ + (double)k*hz;
218  if (k==0 || k==NZ || j==0 || i==0 || j==NY || i==NX)
219  {
220  nodeOnBoundary(inode)=1;
221  }
222  else
223  {
224  nodeOnBoundary(inode)=0;
225  }
226  inode++;
227  }
228  }
229  }
230 #define DUMP_DATA
231 #ifdef DUMP_DATA
232  // Print nodal coords
233  ofstream fcoordout("coords.dat");
234  for (int i=0; i<numNodes; i++) {
235  fcoordout << nodeCoord(i,0) <<" ";
236  fcoordout << nodeCoord(i,1) <<" ";
237  fcoordout << nodeCoord(i,2) <<"\n";
238  }
239  fcoordout.close();
240 #endif
241 
242 
243  // Element to Node map
244  // We'll keep it around, but this is only the DOFMap if you are in the lowest order case.
245  FieldContainer<int> elemToNode(numElems, numNodesPerElem);
246  int ielem = 0;
247  for (int k=0; k<NZ; k++)
248  {
249  for (int j=0; j<NY; j++)
250  {
251  for (int i=0; i<NX; i++)
252  {
253  elemToNode(ielem,0) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
254  elemToNode(ielem,1) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
255  elemToNode(ielem,2) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
256  elemToNode(ielem,3) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
257  elemToNode(ielem,4) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i;
258  elemToNode(ielem,5) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1;
259  elemToNode(ielem,6) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1;
260  elemToNode(ielem,7) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i;
261  ielem++;
262  }
263  }
264  }
265 #ifdef DUMP_DATA
266  // Output connectivity
267  ofstream fe2nout("elem2node.dat");
268  for (int k=0;k<NZ;k++)
269  {
270  for (int j=0; j<NY; j++)
271  {
272  for (int i=0; i<NX; i++)
273  {
274  int ielem = i + j * NX + k * NY * NY;
275  for (int m=0; m<numNodesPerElem; m++)
276  {
277  fe2nout << elemToNode(ielem,m) <<" ";
278  }
279  fe2nout <<"\n";
280  }
281  }
282  }
283  fe2nout.close();
284 #endif
285 
286  // ************************************ CUBATURE **************************************
287  *outStream << "Getting cubature ... \n\n";
288 
289  // Get numerical integration points and weights
290  DefaultCubatureFactory<double> cubFactory;
291  int cubDegree = 2*deg;
292  Teuchos::RCP<Cubature<double> > quadCub = cubFactory.create(hex_8, cubDegree);
293 
294  int cubDim = quadCub->getDimension();
295  int numCubPoints = quadCub->getNumPoints();
296 
297  FieldContainer<double> cubPoints(numCubPoints, cubDim);
298  FieldContainer<double> cubWeights(numCubPoints);
299 
300  quadCub->getCubature(cubPoints, cubWeights);
301 
302 
303  // ************************************** BASIS ***************************************
304 
305  *outStream << "Getting basis ... \n\n";
306 
307  // Define basis
308  Basis_HGRAD_HEX_Cn_FEM<double, FieldContainer<double> > quadHGradBasis(deg,POINTTYPE_SPECTRAL);
309  int numFieldsG = quadHGradBasis.getCardinality();
310  FieldContainer<double> quadGVals(numFieldsG, numCubPoints);
311  FieldContainer<double> quadGrads(numFieldsG, numCubPoints, spaceDim);
312 
313  // Evaluate basis values and gradients at cubature points
314  quadHGradBasis.getValues(quadGVals, cubPoints, OPERATOR_VALUE);
315  quadHGradBasis.getValues(quadGrads, cubPoints, OPERATOR_GRAD);
316 
317  // create the local-global mapping
318  FieldContainer<int> ltgMapping(numElems,numFieldsG);
319  const int numDOF = (NX*deg+1)*(NY*deg+1)*(NZ*deg+1);
320  ielem=0;
321  for (int k=0;k<NZ;k++)
322  {
323  for (int j=0;j<NY;j++)
324  {
325  for (int i=0;i<NX;i++)
326  {
327  const int start = k * ( NY * deg + 1 ) * ( NX * deg + 1 ) + j * ( NX * deg + 1 ) + i * deg;
328  // loop over local dof on this cell
329  int local_dof_cur=0;
330  for (int kloc=0;kloc<=deg;kloc++)
331  {
332  for (int jloc=0;jloc<=deg;jloc++)
333  {
334  for (int iloc=0;iloc<=deg;iloc++)
335  {
336  ltgMapping(ielem,local_dof_cur) = start
337  + kloc * ( NX * deg + 1 ) * ( NY * deg + 1 )
338  + jloc * ( NX * deg + 1 )
339  + iloc;
340  local_dof_cur++;
341  }
342  }
343  }
344  ielem++;
345  }
346  }
347  }
348 #ifdef DUMP_DATA
349  // Output ltg mapping
350  ielem = 0;
351  ofstream ltgout("ltg.dat");
352  for (int k=0;k<NZ;k++)
353  {
354  for (int j=0; j<NY; j++)
355  {
356  for (int i=0; i<NX; i++)
357  {
358  int ielem = i + j * NX + k * NX * NY;
359  for (int m=0; m<numFieldsG; m++)
360  {
361  ltgout << ltgMapping(ielem,m) <<" ";
362  }
363  ltgout <<"\n";
364  }
365  }
366  }
367  ltgout.close();
368 #endif
369 
370  // ********** DECLARE GLOBAL OBJECTS *************
371  Epetra_SerialComm Comm;
372  Epetra_Map globalMapG(numDOF, 0, Comm);
373  Epetra_FEVector u(globalMapG); u.Random();
374  Epetra_FEVector Ku(globalMapG);
375 
376  // time the instantiation
377 // Epetra_Time instantiateTimer(Comm);
378 // Epetra_FECrsMatrix StiffMatrix(Copy,globalMapG,8*numFieldsG);
379 // const double instantiateTime = instantiateTimer.ElapsedTime();
380 
381 
382  // ********** CONSTRUCT AND INSERT LOCAL STIFFNESS MATRICES ***********
383  *outStream << "Building local stiffness matrices...\n\n";
385  typedef FunctionSpaceTools fst;
386  int numCells = numElems;
387 
388  // vertices
389  FieldContainer<double> cellVertices(numCells,numNodesPerElem,spaceDim);
390 
391  // jacobian information
392  FieldContainer<double> cellJacobian(numCells,numCubPoints,spaceDim,spaceDim);
393  FieldContainer<double> cellJacobInv(numCells,numCubPoints,spaceDim,spaceDim);
394  FieldContainer<double> cellJacobDet(numCells,numCubPoints);
395 
396  // element stiffness matrices and supporting storage space
397  FieldContainer<double> localStiffMatrices(numCells, numFieldsG, numFieldsG);
398  FieldContainer<double> transformedBasisGradients(numCells,numFieldsG,numCubPoints,spaceDim);
399  FieldContainer<double> weightedTransformedBasisGradients(numCells,numFieldsG,numCubPoints,spaceDim);
400  FieldContainer<double> weightedMeasure(numCells, numCubPoints);
401 
402 
403  // get vertices of cells (for computing Jacobians)
404  for (int i=0;i<numElems;i++)
405  {
406  for (int j=0;j<numNodesPerElem;j++)
407  {
408  const int nodeCur = elemToNode(i,j);
409  for (int k=0;k<spaceDim;k++)
410  {
411  cellVertices(i,j,k) = nodeCoord(nodeCur,k);
412  }
413  }
414  }
415 
416  Epetra_Time localConstructTimer( Comm );
417 
418  // jacobian evaluation
419  CellTools::setJacobian(cellJacobian,cubPoints,cellVertices,hex_8);
420  CellTools::setJacobianInv(cellJacobInv, cellJacobian );
421  CellTools::setJacobianDet(cellJacobDet, cellJacobian );
422 
423  // transform reference element gradients to each cell
424  fst::HGRADtransformGRAD<double>(transformedBasisGradients, cellJacobInv, quadGrads);
425 
426  // compute weighted measure
427  fst::computeCellMeasure<double>(weightedMeasure, cellJacobDet, cubWeights);
428 
429  // multiply values with weighted measure
430  fst::multiplyMeasure<double>(weightedTransformedBasisGradients,
431  weightedMeasure, transformedBasisGradients);
432 
433  // integrate to compute element stiffness matrix
434  fst::integrate<double>(localStiffMatrices,
435  transformedBasisGradients, weightedTransformedBasisGradients , COMP_BLAS);
436 
437  const double localConstructTime = localConstructTimer.ElapsedTime();
438 
439 
440  Epetra_Time insertionTimer(Comm);
441 
442  vector<map<int,double> > mat(numDOF);
443 
444 
445 
446  // *** Element loop ***
447  for (int el=0; el<numElems; el++)
448  {
449  for (int i=0;i<numFieldsG;i++) // local rows
450  {
451  const int glob_row = ltgMapping(el,i);
452  map<int,double> & cur_row = mat[glob_row];
453 
454  for (int j=0;j<numFieldsG;j++) // local columns
455  {
456  const int glob_col = ltgMapping(el,j);
457  const double cur_val = localStiffMatrices(el,i,j);
458  map<int,double>::iterator it = cur_row.find( glob_col );
459  if (it != cur_row.end()) // current column already in row
460  {
461  it->second += cur_val;
462  }
463  else
464  {
465  cur_row[glob_col] = cur_val;
466  }
467  }
468  }
469  }
470  //StiffMatrix.GlobalAssemble(); StiffMatrix.FillComplete();
471  const double insertionTime = insertionTimer.ElapsedTime( );
472 
473 // *outStream << "Time to instantiate global stiffness matrix: " << instantiateTime << "\n";
474  *outStream << "Time to build local matrices (including Jacobian computation): "<< localConstructTime << "\n";
475  *outStream << "Time to assemble global matrix from local matrices: " << insertionTime << "\n";
476  *outStream << "Total construction time: " << localConstructTime + insertionTime << "\n";
477 
478 // Epetra_Time applyTimer(Comm);
479 // StiffMatrix.Apply(u,Ku);
480 // const double multTime = applyTimer.ElapsedTime();
481 // *outStream << "Time to multiply onto a vector: " << multTime << "\n";
482 
483  *outStream << "End Result: TEST PASSED\n";
484 
485  // reset format state of std::cout
486  std::cout.copyfmt(oldFormatState);
487 
488  return 0;
489 }
490 
Header file for the Intrepid::CellTools class.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_HEX_Cn_FEM class.
Defines expert-level interfaces for the evaluation of functions and operators in physical space (supp...
Intrepid utilities.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell...
Header file for the Intrepid::FunctionSpaceTools class.
A factory class that generates specific instances of cubatures.
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > &degree)
Factory method.
A stateless class for operations on cell data. Provides methods for: