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