Intrepid
example_07.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"
90 #include "Intrepid_Utils.hpp"
91 
92 // Epetra includes
93 #include "Epetra_Time.h"
94 #include "Epetra_Map.h"
95 #include "Epetra_FECrsMatrix.h"
96 #include "Epetra_FEVector.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 
104 // Shards includes
105 #include "Shards_CellTopology.hpp"
106 
107 // EpetraExt includes
108 #include "EpetraExt_RowMatrixOut.h"
109 #include "EpetraExt_MultiVectorOut.h"
110 
111 using namespace std;
112 using namespace Intrepid;
113 
114 // Functions to evaluate exact solution and derivatives
115 double evalu(double & x, double & y, double & z);
116 int evalGradu(double & x, double & y, double & z, double & gradu1, double & gradu2, double & gradu3);
117 double evalDivGradu(double & x, double & y, double & z);
118 
119 int main(int argc, char *argv[]) {
120 
121  //Check number of arguments
122  if (argc < 4) {
123  std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
124  std::cout <<"Usage:\n\n";
125  std::cout <<" ./Intrepid_example_Drivers_Example_05.exe deg NX NY verbose\n\n";
126  std::cout <<" where \n";
127  std::cout <<" int deg - polynomial degree to be used (assumed > 1) \n";
128  std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n";
129  std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n";
130  std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
131  exit(1);
132  }
133 
134  // This little trick lets us print to std::cout only if
135  // a (dummy) command-line argument is provided.
136  int iprint = argc - 1;
137  Teuchos::RCP<std::ostream> outStream;
138  Teuchos::oblackholestream bhs; // outputs nothing
139  if (iprint > 2)
140  outStream = Teuchos::rcp(&std::cout, false);
141  else
142  outStream = Teuchos::rcp(&bhs, false);
143 
144  // Save the format state of the original std::cout.
145  Teuchos::oblackholestream oldFormatState;
146  oldFormatState.copyfmt(std::cout);
147 
148  *outStream \
149  << "===============================================================================\n" \
150  << "| |\n" \
151  << "| Example: Generate Stiffness Matrix and Right Hand Side Vector for |\n" \
152  << "| Poisson Equation on Quadrilateral Mesh |\n" \
153  << "| |\n" \
154  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
155  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
156  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
157  << "| |\n" \
158  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
159  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
160  << "| |\n" \
161  << "===============================================================================\n";
162 
163 
164  // ************************************ GET INPUTS **************************************
165 
166  int deg = atoi(argv[1]); // polynomial degree to use
167  int NX = atoi(argv[2]); // num intervals in x direction (assumed box domain, 0,1)
168  int NY = atoi(argv[3]); // 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 quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
176 
177  // Get dimensions
178  int numNodesPerElem = quad_4.getNodeCount();
179  int spaceDim = quad_4.getDimension();
180 
181  // *********************************** GENERATE MESH ************************************
182 
183  *outStream << "Generating mesh ... \n\n";
184 
185  *outStream << " NX" << " NY\n";
186  *outStream << std::setw(5) << NX <<
187  std::setw(5) << NY << "\n\n";
188 
189  // Print mesh information
190  int numElems = NX*NY;
191  int numNodes = (NX+1)*(NY+1);
192  *outStream << " Number of Elements: " << numElems << " \n";
193  *outStream << " Number of Nodes: " << numNodes << " \n\n";
194 
195  // Square
196  double leftX = 0.0, rightX = 1.0;
197  double leftY = 0.0, rightY = 1.0;
198 
199  // Mesh spacing
200  double hx = (rightX-leftX)/((double)NX);
201  double hy = (rightY-leftY)/((double)NY);
202 
203  // Get nodal coordinates
204  FieldContainer<double> nodeCoord(numNodes, spaceDim);
205  FieldContainer<int> nodeOnBoundary(numNodes);
206  int inode = 0;
207  for (int j=0; j<NY+1; j++) {
208  for (int i=0; i<NX+1; i++) {
209  nodeCoord(inode,0) = leftX + (double)i*hx;
210  nodeCoord(inode,1) = leftY + (double)j*hy;
211  if (j==0 || i==0 || j==NY || i==NX){
212  nodeOnBoundary(inode)=1;
213  }
214  else {
215  nodeOnBoundary(inode)=0;
216  }
217  inode++;
218  }
219  }
220 #define DUMP_DATA
221 #ifdef DUMP_DATA
222  // Print nodal coords
223  ofstream fcoordout("coords.dat");
224  for (int i=0; i<numNodes; i++) {
225  fcoordout << nodeCoord(i,0) <<" ";
226  fcoordout << nodeCoord(i,1) <<"\n";
227  }
228  fcoordout.close();
229 #endif
230 
231 
232  // Element to Node map
233  // We'll keep it around, but this is only the DOFMap if you are in the lowest order case.
234  FieldContainer<int> elemToNode(numElems, numNodesPerElem);
235  int ielem = 0;
236  for (int j=0; j<NY; j++) {
237  for (int i=0; i<NX; i++) {
238  elemToNode(ielem,0) = (NX + 1)*j + i;
239  elemToNode(ielem,1) = (NX + 1)*j + i + 1;
240  elemToNode(ielem,2) = (NX + 1)*(j + 1) + i + 1;
241  elemToNode(ielem,3) = (NX + 1)*(j + 1) + i;
242  ielem++;
243  }
244  }
245 #ifdef DUMP_DATA
246  // Output connectivity
247  ofstream fe2nout("elem2node.dat");
248  for (int j=0; j<NY; j++) {
249  for (int i=0; i<NX; i++) {
250  int ielem = i + j * NX;
251  for (int m=0; m<numNodesPerElem; m++){
252  fe2nout << elemToNode(ielem,m) <<" ";
253  }
254  fe2nout <<"\n";
255  }
256  }
257  fe2nout.close();
258 #endif
259 
260 
261  // ************************************ CUBATURE **************************************
262  *outStream << "Getting cubature ... \n\n";
263 
264  // Get numerical integration points and weights
265  DefaultCubatureFactory<double> cubFactory;
266  int cubDegree = 2*deg;
267  Teuchos::RCP<Cubature<double> > quadCub = cubFactory.create(quad_4, cubDegree);
268 
269  int cubDim = quadCub->getDimension();
270  int numCubPoints = quadCub->getNumPoints();
271 
272  FieldContainer<double> cubPoints(numCubPoints, cubDim);
273  FieldContainer<double> cubWeights(numCubPoints);
274 
275  quadCub->getCubature(cubPoints, cubWeights);
276 
277 
278  // ************************************** BASIS ***************************************
279 
280  *outStream << "Getting basis ... \n\n";
281 
282  // Define basis
283  Basis_HGRAD_QUAD_Cn_FEM<double, FieldContainer<double> > quadHGradBasis(deg,POINTTYPE_SPECTRAL);
284  int numFieldsG = quadHGradBasis.getCardinality();
285  FieldContainer<double> quadGVals(numFieldsG, numCubPoints);
286  FieldContainer<double> quadGrads(numFieldsG, numCubPoints, spaceDim);
287 
288  // Evaluate basis values and gradients at cubature points
289  quadHGradBasis.getValues(quadGVals, cubPoints, OPERATOR_VALUE);
290  quadHGradBasis.getValues(quadGrads, cubPoints, OPERATOR_GRAD);
291 
292  // create the local-global mapping for higher order elements
293  FieldContainer<int> ltgMapping(numElems,numFieldsG);
294  const int numDOF = (NX*deg+1)*(NY*deg+1);
295  ielem=0;
296  for (int j=0;j<NY;j++) {
297  for (int i=0;i<NX;i++) {
298  const int start = deg * j * ( NX * deg + 1 ) + i * deg;
299  // loop over local dof on this cell
300  int local_dof_cur=0;
301  for (int vertical=0;vertical<=deg;vertical++) {
302  for (int horizontal=0;horizontal<=deg;horizontal++) {
303  ltgMapping(ielem,local_dof_cur) = start + vertical*(NX*deg+1)+horizontal;
304  local_dof_cur++;
305  }
306  }
307  ielem++;
308  }
309  }
310 #ifdef DUMP_DATA
311  // Output ltg mapping
312  ofstream ltgout("ltg.dat");
313  for (int j=0; j<NY; j++) {
314  for (int i=0; i<NX; i++) {
315  int ielem = i + j * NX;
316  for (int m=0; m<numFieldsG; m++){
317  ltgout << ltgMapping(ielem,m) <<" ";
318  }
319  ltgout <<"\n";
320  }
321  }
322  ltgout.close();
323 #endif
324 
325  // ******** CREATE A SINGLE STIFFNESS MATRIX, WHICH IS REPLICATED ON ALL ELEMENTS *********
326  *outStream << "Building stiffness matrix and right hand side ... \n\n";
327 
328  // Settings and data structures for mass and stiffness matrices
330  typedef FunctionSpaceTools fst;
331  int numCells = 1;
332 
333  // Container for nodes
334  FieldContainer<double> refQuadNodes(numCells, numNodesPerElem, spaceDim);
335  // Containers for Jacobian
336  FieldContainer<double> refQuadJacobian(numCells, numCubPoints, spaceDim, spaceDim);
337  FieldContainer<double> refQuadJacobInv(numCells, numCubPoints, spaceDim, spaceDim);
338  FieldContainer<double> refQuadJacobDet(numCells, numCubPoints);
339  // Containers for element HGRAD stiffness matrix
340  FieldContainer<double> localStiffMatrix(numCells, numFieldsG, numFieldsG);
341  FieldContainer<double> weightedMeasure(numCells, numCubPoints);
342  FieldContainer<double> quadGradsTransformed(numCells, numFieldsG, numCubPoints, spaceDim);
343  FieldContainer<double> quadGradsTransformedWeighted(numCells, numFieldsG, numCubPoints, spaceDim);
344  // Containers for right hand side vectors
345  FieldContainer<double> rhsData(numCells, numCubPoints);
346  FieldContainer<double> localRHS(numCells, numFieldsG);
347  FieldContainer<double> quadGValsTransformed(numCells, numFieldsG, numCubPoints);
348  FieldContainer<double> quadGValsTransformedWeighted(numCells, numFieldsG, numCubPoints);
349  // Container for cubature points in physical space
350  FieldContainer<double> physCubPoints(numCells, numCubPoints, cubDim);
351 
352  // Global arrays in Epetra format
353  // we will explicitly build the sparsity pattern before instantiating the matrix later.
354  Epetra_SerialComm Comm;
355  Epetra_Map globalMapG(numDOF, 0, Comm);
356  Epetra_FEVector u(globalMapG);
357  Epetra_FEVector Ku(globalMapG);
358  u.Random();
359 
360  // ************************** Compute element HGrad stiffness matrices *******************************
361  refQuadNodes(0,0,0) = 0.0;
362  refQuadNodes(0,0,1) = 0.0;
363  refQuadNodes(0,1,0) = hx;
364  refQuadNodes(0,1,1) = 0.0;
365  refQuadNodes(0,2,0) = hx;
366  refQuadNodes(0,2,1) = hy;
367  refQuadNodes(0,3,0) = 0.0;
368  refQuadNodes(0,3,1) = hy;
369 
370  // Compute cell Jacobians, their inverses and their determinants
371  CellTools::setJacobian(refQuadJacobian, cubPoints, refQuadNodes, quad_4);
372  CellTools::setJacobianInv(refQuadJacobInv, refQuadJacobian );
373  CellTools::setJacobianDet(refQuadJacobDet, refQuadJacobian );
374 
375  // transform from [-1,1]^2 to [0,hx]x[0,hy]
376  fst::HGRADtransformGRAD<double>(quadGradsTransformed, refQuadJacobInv, quadGrads);
377 
378  // compute weighted measure
379  fst::computeCellMeasure<double>(weightedMeasure, refQuadJacobDet, cubWeights);
380 
381  // multiply values with weighted measure
382  fst::multiplyMeasure<double>(quadGradsTransformedWeighted,
383  weightedMeasure, quadGradsTransformed);
384 
385  // integrate to compute element stiffness matrix
386  fst::integrate<double>(localStiffMatrix,
387  quadGradsTransformed, quadGradsTransformedWeighted, COMP_BLAS);
388 
389  Epetra_Time graphTimer(Comm);
390  Epetra_CrsGraph grph( Copy , globalMapG , 4 * numFieldsG );
391  for (int k=0;k<numElems;k++)
392  {
393  for (int i=0;i<numFieldsG;i++)
394  {
395  grph.InsertGlobalIndices(ltgMapping(k,i),numFieldsG,&ltgMapping(k,0));
396  }
397  }
398  grph.FillComplete();
399  const double graphTime = graphTimer.ElapsedTime();
400  std::cout << "Graph computed in " << graphTime << "\n";
401 
402  Epetra_Time instantiateTimer( Comm );
403  Epetra_FECrsMatrix StiffMatrix( Copy , grph );
404  const double instantiateTime = instantiateTimer.ElapsedTime( );
405  std::cout << "Matrix instantiated in " << instantiateTime << "\n";
406 
407  Epetra_Time assemblyTimer(Comm);
408 
409  // *** Element loop ***
410  for (int k=0; k<numElems; k++)
411  {
412  // assemble into global matrix
413  StiffMatrix.InsertGlobalValues(numFieldsG,&ltgMapping(k,0),numFieldsG,&ltgMapping(k,0),&localStiffMatrix(0,0,0));
414 
415  }
416 
417 
418  // Assemble global matrices
419  StiffMatrix.GlobalAssemble(); StiffMatrix.FillComplete();
420 
421  double assembleTime = assemblyTimer.ElapsedTime();
422  std::cout << "Time to insert reference element matrix into global matrix: " << assembleTime << std::endl;
423  std::cout << "Total matrix construction time: " << assembleTime + instantiateTime + graphTime << "\n";
424  std::cout << "There are " << StiffMatrix.NumGlobalNonzeros() << " nonzeros in the matrix.\n";
425  std::cout << "There are " << numDOF << " global degrees of freedom.\n";
426 
427  Epetra_Time multTimer(Comm);
428  StiffMatrix.Apply(u,Ku);
429  double multTime = multTimer.ElapsedTime();
430  std::cout << "Time to apply: " << multTime << std::endl;
431 
432 
433 #ifdef DUMP_DATA
434  // Dump matrices to disk
435 // EpetraExt::RowMatrixToMatlabFile("stiff_matrix.dat",StiffMatrix);
436 // EpetraExt::MultiVectorToMatrixMarketFile("rhs_vector.dat",rhs,0,0,false);
437 #endif
438 
439  std::cout << "End Result: TEST PASSED\n";
440 
441  // reset format state of std::cout
442  std::cout.copyfmt(oldFormatState);
443 
444  return 0;
445 }
446 
447 
448 // Calculates value of exact solution u
449  double evalu(double & x, double & y, double & z)
450  {
451  /*
452  // function1
453  double exactu = sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
454  */
455 
456  // function2
457  double exactu = sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z);
458 
459  return exactu;
460  }
461 
462 // Calculates gradient of exact solution u
463  int evalGradu(double & x, double & y, double & z, double & gradu1, double & gradu2, double & gradu3)
464  {
465  /*
466  // function 1
467  gradu1 = M_PI*cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
468  gradu2 = M_PI*sin(M_PI*x)*cos(M_PI*y)*sin(M_PI*z);
469  gradu3 = M_PI*sin(M_PI*x)*sin(M_PI*y)*cos(M_PI*z);
470  */
471 
472  // function2
473  gradu1 = (M_PI*cos(M_PI*x)+sin(M_PI*x))
474  *sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z);
475  gradu2 = (M_PI*cos(M_PI*y)+sin(M_PI*y))
476  *sin(M_PI*x)*sin(M_PI*z)*exp(x+y+z);
477  gradu3 = (M_PI*cos(M_PI*z)+sin(M_PI*z))
478  *sin(M_PI*x)*sin(M_PI*y)*exp(x+y+z);
479 
480  return 0;
481  }
482 
483 // Calculates Laplacian of exact solution u
484  double evalDivGradu(double & x, double & y, double & z)
485  {
486  /*
487  // function 1
488  double divGradu = -3.0*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
489  */
490 
491  // function 2
492  double divGradu = -3.0*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z)
493  + 2.0*M_PI*cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z)
494  + 2.0*M_PI*cos(M_PI*y)*sin(M_PI*x)*sin(M_PI*z)*exp(x+y+z)
495  + 2.0*M_PI*cos(M_PI*z)*sin(M_PI*x)*sin(M_PI*y)*exp(x+y+z)
496  + 3.0*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z);
497 
498  return divGradu;
499  }
500 
Header file for the Intrepid::CellTools class.
Header file for utility class to provide multidimensional containers.
Header file for utility class to provide array tools, such as tensor contractions, etc.
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.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
Header file for the Intrepid::HGRAD_QUAD_Cn_FEM 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: