Intrepid
example_08.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 
81 // Intrepid includes
84 #include "Intrepid_CellTools.hpp"
85 #include "Intrepid_ArrayTools.hpp"
89 #include "Intrepid_Utils.hpp"
90 
91 // Epetra includes
92 #include "Epetra_Time.h"
93 #include "Epetra_Map.h"
94 #include "Epetra_FECrsMatrix.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 
103 // Shards includes
104 #include "Shards_CellTopology.hpp"
105 
106 // EpetraExt includes
107 #include "EpetraExt_RowMatrixOut.h"
108 #include "EpetraExt_MultiVectorOut.h"
109 
110 using namespace std;
111 using namespace Intrepid;
112 
113 // Functions to evaluate exact solution and derivatives
114 double evalu(double & x, double & y, double & z);
115 int evalGradu(double & x, double & y, double & z, double & gradu1, double & gradu2, double & gradu3);
116 double evalDivGradu(double & x, double & y, double & z);
117 
118 int main(int argc, char *argv[]) {
119 
120  //Check number of arguments
121  if (argc < 4) {
122  std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
123  std::cout <<"Usage:\n\n";
124  std::cout <<" ./Intrepid_example_Drivers_Example_05.exe deg NX NY verbose\n\n";
125  std::cout <<" where \n";
126  std::cout <<" int deg - polynomial degree to be used (assumed > 1) \n";
127  std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n";
128  std::cout <<" int NY - 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: Generate Stiffness Matrix and Right Hand Side Vector for |\n" \
151  << "| Poisson Equation on Quadrilateral 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 
169 
170  // *********************************** CELL TOPOLOGY **********************************
171 
172  // Get cell topology for base hexahedron
173  typedef shards::CellTopology CellTopology;
174  CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
175 
176  // Get dimensions
177  int numNodesPerElem = quad_4.getNodeCount();
178  int spaceDim = quad_4.getDimension();
179 
180  // *********************************** GENERATE MESH ************************************
181 
182  *outStream << "Generating mesh ... \n\n";
183 
184  *outStream << " NX" << " NY\n";
185  *outStream << std::setw(5) << NX <<
186  std::setw(5) << NY << "\n\n";
187 
188  // Print mesh information
189  int numElems = NX*NY;
190  int numNodes = (NX+1)*(NY+1);
191  *outStream << " Number of Elements: " << numElems << " \n";
192  *outStream << " Number of Nodes: " << numNodes << " \n\n";
193 
194  // Square
195  double leftX = 0.0, rightX = 1.0;
196  double leftY = 0.0, rightY = 1.0;
197 
198  // Mesh spacing
199  double hx = (rightX-leftX)/((double)NX);
200  double hy = (rightY-leftY)/((double)NY);
201 
202  // Get nodal coordinates
203  FieldContainer<double> nodeCoord(numNodes, spaceDim);
204  FieldContainer<int> nodeOnBoundary(numNodes);
205  int inode = 0;
206  for (int j=0; j<NY+1; j++) {
207  for (int i=0; i<NX+1; i++) {
208  nodeCoord(inode,0) = leftX + (double)i*hx;
209  nodeCoord(inode,1) = leftY + (double)j*hy;
210  if (j==0 || i==0 || j==NY || i==NX){
211  nodeOnBoundary(inode)=1;
212  }
213  else {
214  nodeOnBoundary(inode)=0;
215  }
216  inode++;
217  }
218  }
219 #define DUMP_DATA
220 #ifdef DUMP_DATA
221  // Print nodal coords
222  ofstream fcoordout("coords.dat");
223  for (int i=0; i<numNodes; i++) {
224  fcoordout << nodeCoord(i,0) <<" ";
225  fcoordout << nodeCoord(i,1) <<"\n";
226  }
227  fcoordout.close();
228 #endif
229 
230 
231  // Element to Node map
232  // We'll keep it around, but this is only the DOFMap if you are in the lowest order case.
233  FieldContainer<int> elemToNode(numElems, numNodesPerElem);
234  int ielem=0;
235  for (int j=0; j<NY; j++) {
236  for (int i=0; i<NX; i++) {
237  elemToNode(ielem,0) = (NX + 1)*j + i;
238  elemToNode(ielem,1) = (NX + 1)*j + i + 1;
239  elemToNode(ielem,2) = (NX + 1)*(j + 1) + i + 1;
240  elemToNode(ielem,3) = (NX + 1)*(j + 1) + i;
241  ielem++;
242  }
243  }
244 #ifdef DUMP_DATA
245  // Output connectivity
246  ofstream fe2nout("elem2node.dat");
247  for (int j=0; j<NY; j++) {
248  for (int i=0; i<NX; i++) {
249  ielem = i + j * NX;
250  for (int m=0; m<numNodesPerElem; m++){
251  fe2nout << elemToNode(ielem,m) <<" ";
252  }
253  fe2nout <<"\n";
254  }
255  }
256  fe2nout.close();
257 #endif
258 
259 
260  // ************************************ CUBATURE **************************************
261  *outStream << "Getting cubature ... \n\n";
262 
263  // Get numerical integration points and weights
264  DefaultCubatureFactory<double> cubFactory;
265  int cubDegree = 2*deg;
266  Teuchos::RCP<Cubature<double> > quadCub = cubFactory.create(quad_4, cubDegree);
267 
268  int cubDim = quadCub->getDimension();
269  int numCubPoints = quadCub->getNumPoints();
270 
271  FieldContainer<double> cubPoints(numCubPoints, cubDim);
272  FieldContainer<double> cubWeights(numCubPoints);
273 
274  quadCub->getCubature(cubPoints, cubWeights);
275 
276 
277  // ************************************** BASIS ***************************************
278 
279  *outStream << "Getting basis ... \n\n";
280 
281  // Define basis
282  Basis_HGRAD_QUAD_Cn_FEM<double, FieldContainer<double> > quadHGradBasis(deg,POINTTYPE_SPECTRAL);
283  int numFieldsG = quadHGradBasis.getCardinality();
284  FieldContainer<double> quadGVals(numFieldsG, numCubPoints);
285  FieldContainer<double> quadGrads(numFieldsG, numCubPoints, spaceDim);
286 
287  // Evaluate basis values and gradients at cubature points
288  quadHGradBasis.getValues(quadGVals, cubPoints, OPERATOR_VALUE);
289  quadHGradBasis.getValues(quadGrads, cubPoints, OPERATOR_GRAD);
290 
291  // create the local-global mapping for higher order elements
292  FieldContainer<int> ltgMapping(numElems,numFieldsG);
293  const int numDOF = (NX*deg+1)*(NY*deg+1);
294  ielem=0;
295  for (int j=0;j<NY;j++) {
296  for (int i=0;i<NX;i++) {
297  const int start = deg * j * ( NX * deg + 1 ) + i * deg;
298  // loop over local dof on this cell
299  int local_dof_cur=0;
300  for (int vertical=0;vertical<=deg;vertical++) {
301  for (int horizontal=0;horizontal<=deg;horizontal++) {
302  ltgMapping(ielem,local_dof_cur) = start + vertical*(NX*deg+1)+horizontal;
303  local_dof_cur++;
304  }
305  }
306  ielem++;
307  }
308  }
309 #ifdef DUMP_DATA
310  // Output ltg mapping
311  ofstream ltgout("ltg.dat");
312  for (int j=0; j<NY; j++) {
313  for (int i=0; i<NX; i++) {
314  ielem = i + j * NX;
315  for (int m=0; m<numFieldsG; m++){
316  ltgout << ltgMapping(ielem,m) <<" ";
317  }
318  ltgout <<"\n";
319  }
320  }
321  ltgout.close();
322 #endif
323 
324  // ******** CREATE ALL LOCAL STIFFNESS MATRICES *********
325  *outStream << "Building stiffness matrix and right hand side ... \n\n";
326 
327  // Settings and data structures for mass and stiffness matrices
329  typedef FunctionSpaceTools fst;
330  int numCells = numElems;
331 
332  // Container for nodes
333  FieldContainer<double> cellVertices(numCells,numNodesPerElem,spaceDim);
334 
335  // Containers for Jacobian
336  FieldContainer<double> cellJacobian(numCells, numCubPoints, spaceDim, spaceDim);
337  FieldContainer<double> cellJacobInv(numCells, numCubPoints, spaceDim, spaceDim);
338  FieldContainer<double> cellJacobDet(numCells, numCubPoints);
339 
340  // Containers for element HGRAD stiffness matrices
341  FieldContainer<double> localStiffMatrices(numCells, numFieldsG, numFieldsG);
342  FieldContainer<double> transformedBasisGradients(numCells,numFieldsG,numCubPoints,spaceDim);
343  FieldContainer<double> weightedTransformedBasisGradients(numCells,numFieldsG,numCubPoints,spaceDim);
344  FieldContainer<double> weightedMeasure(numCells, numCubPoints);
345 
346 
347  // Global arrays in Epetra format
348  Epetra_SerialComm Comm;
349  Epetra_Map globalMapG(numDOF, 0, Comm);
350 
351  Epetra_Time graphTimer(Comm);
352  Epetra_CrsGraph grph( Copy , globalMapG , 4 * numFieldsG );
353  for (int k=0;k<numElems;k++)
354  {
355  for (int i=0;i<numFieldsG;i++)
356  {
357  grph.InsertGlobalIndices(ltgMapping(k,i),numFieldsG,&ltgMapping(k,0));
358  }
359  }
360  grph.FillComplete();
361 
362  const double graphTime = graphTimer.ElapsedTime();
363  std::cout << "Graph computed in " << graphTime << "\n";
364 
365  Epetra_Time instantiateTimer( Comm );
366  Epetra_FECrsMatrix StiffMatrix( Copy , grph );
367  const double instantiateTime = instantiateTimer.ElapsedTime( );
368  std::cout << "Matrix instantiated in " << instantiateTime << "\n";
369 
370  Epetra_FEVector u(globalMapG);
371  Epetra_FEVector Ku(globalMapG);
372 
373  u.Random();
374 
375  // ************************** Compute element HGrad stiffness matrices *******************************
376  // Get vertices of all the cells
377 
378  for (int i=0;i<numElems;i++)
379  {
380  for (int j=0;j<4;j++)
381  {
382  const int nodeCur = elemToNode(i,j);
383  for (int k=0;k<spaceDim;k++)
384  {
385  cellVertices(i,j,k) = nodeCoord(nodeCur,k);
386  }
387  }
388  }
389 
390  Epetra_Time localConstructTimer(Comm);
391 
392  // Compute all cell Jacobians, their inverses and their determinants
393  CellTools::setJacobian(cellJacobian, cubPoints, cellVertices, quad_4);
394  CellTools::setJacobianInv(cellJacobInv, cellJacobian );
395  CellTools::setJacobianDet(cellJacobDet, cellJacobian );
396 
397  // transform reference element gradients to each cell
398  fst::HGRADtransformGRAD<double>(transformedBasisGradients, cellJacobInv, quadGrads);
399 
400  // compute weighted measure
401  fst::computeCellMeasure<double>(weightedMeasure, cellJacobDet, cubWeights);
402 
403  // multiply values with weighted measure
404  fst::multiplyMeasure<double>(weightedTransformedBasisGradients,
405  weightedMeasure, transformedBasisGradients);
406 
407  // integrate to compute element stiffness matrix
408  fst::integrate<double>(localStiffMatrices,
409  transformedBasisGradients, weightedTransformedBasisGradients , COMP_BLAS);
410 
411  const double localConstructTime = localConstructTimer.ElapsedTime();
412  std::cout << "Time to build local matrices (including Jacobian computation): "<< localConstructTime << "\n";
413 
414  Epetra_Time insertionTimer(Comm);
415 
416  // *** Element loop ***
417  for (int k=0; k<numElems; k++)
418  {
419  // assemble into global matrix
420  StiffMatrix.InsertGlobalValues(numFieldsG,&ltgMapping(k,0),numFieldsG,&ltgMapping(k,0),&localStiffMatrices(k,0,0));
421 
422  }
423  StiffMatrix.GlobalAssemble(); StiffMatrix.FillComplete();
424  const double insertionTime = insertionTimer.ElapsedTime( );
425 
426  std::cout << "Time to assemble global matrix from local matrices: " << insertionTime << "\n";
427 
428 
429 #ifdef DUMP_DATA
430  // Dump matrices to disk
431 // EpetraExt::RowMatrixToMatlabFile("stiff_matrix.dat",StiffMatrix);
432 // EpetraExt::MultiVectorToMatrixMarketFile("rhs_vector.dat",rhs,0,0,false);
433 #endif
434 
435 
436  std::cout << "End Result: TEST PASSED\n";
437 
438  // reset format state of std::cout
439  std::cout.copyfmt(oldFormatState);
440 
441  return 0;
442 }
443 
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: