Intrepid
example_06.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_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_06.exe deg NX NY 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 <<" verbose (optional) - any character, indicates verbose output \n\n";
125  exit(1);
126  }
127 
128  // This little trick lets us print to std::cout only if
129  // a (dummy) command-line argument is provided.
130  int iprint = argc - 1;
131  Teuchos::RCP<std::ostream> outStream;
132  Teuchos::oblackholestream bhs; // outputs nothing
133  if (iprint > 2)
134  outStream = Teuchos::rcp(&std::cout, false);
135  else
136  outStream = Teuchos::rcp(&bhs, false);
137 
138  // Save the format state of the original std::cout.
139  Teuchos::oblackholestream oldFormatState;
140  oldFormatState.copyfmt(std::cout);
141 
142  *outStream \
143  << "===============================================================================\n" \
144  << "| |\n" \
145  << "| Example: Apply Stiffness Matrix for |\n" \
146  << "| Poisson Equation on Quadrilateral Mesh |\n" \
147  << "| |\n" \
148  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
149  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
150  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
151  << "| |\n" \
152  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
153  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
154  << "| |\n" \
155  << "===============================================================================\n";
156 
157 
158  // ************************************ GET INPUTS **************************************
159 
160  int deg = atoi(argv[1]); // polynomial degree to use
161  int NX = atoi(argv[2]); // num intervals in x direction (assumed box domain, 0,1)
162  int NY = atoi(argv[3]); // num intervals in y direction (assumed box domain, 0,1)
163 
164 
165  // *********************************** CELL TOPOLOGY **********************************
166 
167  // Get cell topology for base hexahedron
168  typedef shards::CellTopology CellTopology;
169  CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
170 
171  // Get dimensions
172  int numNodesPerElem = quad_4.getNodeCount();
173  int spaceDim = quad_4.getDimension();
174 
175  // *********************************** GENERATE MESH ************************************
176 
177  *outStream << "Generating mesh ... \n\n";
178 
179  *outStream << " NX" << " NY\n";
180  *outStream << std::setw(5) << NX <<
181  std::setw(5) << NY << "\n\n";
182 
183  // Print mesh information
184  int numElems = NX*NY;
185  int numNodes = (NX+1)*(NY+1);
186  *outStream << " Number of Elements: " << numElems << " \n";
187  *outStream << " Number of Nodes: " << numNodes << " \n\n";
188 
189  // Square
190  double leftX = 0.0, rightX = 1.0;
191  double leftY = 0.0, rightY = 1.0;
192 
193  // Mesh spacing
194  double hx = (rightX-leftX)/((double)NX);
195  double hy = (rightY-leftY)/((double)NY);
196 
197  // Get nodal coordinates
198  FieldContainer<double> nodeCoord(numNodes, spaceDim);
199  FieldContainer<int> nodeOnBoundary(numNodes);
200  int inode = 0;
201  for (int j=0; j<NY+1; j++) {
202  for (int i=0; i<NX+1; i++) {
203  nodeCoord(inode,0) = leftX + (double)i*hx;
204  nodeCoord(inode,1) = leftY + (double)j*hy;
205  if (j==0 || i==0 || j==NY || i==NX){
206  nodeOnBoundary(inode)=1;
207  }
208  else {
209  nodeOnBoundary(inode)=0;
210  }
211  inode++;
212  }
213  }
214 #define DUMP_DATA
215 #ifdef DUMP_DATA
216  // Print nodal coords
217  ofstream fcoordout("coords.dat");
218  for (int i=0; i<numNodes; i++) {
219  fcoordout << nodeCoord(i,0) <<" ";
220  fcoordout << nodeCoord(i,1) <<"\n";
221  }
222  fcoordout.close();
223 #endif
224 
225 
226  // Element to Node map
227  // We'll keep it around, but this is only the DOFMap if you are in the lowest order case.
228  FieldContainer<int> elemToNode(numElems, numNodesPerElem);
229  int ielem = 0;
230  for (int j=0; j<NY; j++) {
231  for (int i=0; i<NX; i++) {
232  elemToNode(ielem,0) = (NX + 1)*j + i;
233  elemToNode(ielem,1) = (NX + 1)*j + i + 1;
234  elemToNode(ielem,2) = (NX + 1)*(j + 1) + i + 1;
235  elemToNode(ielem,3) = (NX + 1)*(j + 1) + i;
236  ielem++;
237  }
238  }
239 #ifdef DUMP_DATA
240  // Output connectivity
241  ofstream fe2nout("elem2node.dat");
242  for (int j=0; j<NY; j++) {
243  for (int i=0; i<NX; i++) {
244  int ielem = i + j * NX;
245  for (int m=0; m<numNodesPerElem; m++){
246  fe2nout << elemToNode(ielem,m) <<" ";
247  }
248  fe2nout <<"\n";
249  }
250  }
251  fe2nout.close();
252 #endif
253 
254  // ************************************ CUBATURE **************************************
255  *outStream << "Getting cubature ... \n\n";
256 
257  // Get numerical integration points and weights
258  DefaultCubatureFactory<double> cubFactory;
259  int cubDegree = 2*deg;
260  Teuchos::RCP<Cubature<double> > quadCub = cubFactory.create(quad_4, cubDegree);
261 
262  int cubDim = quadCub->getDimension();
263  int numCubPoints = quadCub->getNumPoints();
264 
265  FieldContainer<double> cubPoints(numCubPoints, cubDim);
266  FieldContainer<double> cubWeights(numCubPoints);
267 
268  quadCub->getCubature(cubPoints, cubWeights);
269 
270 
271  // ************************************** BASIS ***************************************
272 
273  *outStream << "Getting basis ... \n\n";
274 
275  // Define basis
276  Basis_HGRAD_QUAD_Cn_FEM<double, FieldContainer<double> > quadHGradBasis(deg,POINTTYPE_SPECTRAL);
277  int numFieldsG = quadHGradBasis.getCardinality();
278  FieldContainer<double> quadGVals(numFieldsG, numCubPoints);
279  FieldContainer<double> quadGrads(numFieldsG, numCubPoints, spaceDim);
280 
281  // Evaluate basis values and gradients at cubature points
282  quadHGradBasis.getValues(quadGVals, cubPoints, OPERATOR_VALUE);
283  quadHGradBasis.getValues(quadGrads, cubPoints, OPERATOR_GRAD);
284 
285  // create the local-global mapping for higher order elements
286  FieldContainer<int> ltgMapping(numElems,numFieldsG);
287  const int numDOF = (NX*deg+1)*(NY*deg+1);
288  ielem=0;
289  for (int j=0;j<NY;j++) {
290  for (int i=0;i<NX;i++) {
291  const int start = deg * j * ( NX * deg + 1 ) + i * deg;
292  // loop over local dof on this cell
293  int local_dof_cur=0;
294  for (int vertical=0;vertical<=deg;vertical++) {
295  for (int horizontal=0;horizontal<=deg;horizontal++) {
296  ltgMapping(ielem,local_dof_cur) = start + vertical*(NX*deg+1)+horizontal;
297  local_dof_cur++;
298  }
299  }
300  ielem++;
301  }
302  }
303 #ifdef DUMP_DATA
304  // Output ltg mapping
305 // ofstream ltgout("ltg.dat");
306 // for (int j=0; j<NY; j++) {
307 // for (int i=0; i<NX; i++) {
308 // int ielem = i + j * NX;
309 // for (int m=0; m<numFieldsG; m++){
310 // ltgout << ltgMapping(ielem,m) <<" ";
311 // }
312 // ltgout <<"\n";
313 // }
314 // }
315 // ltgout.close();
316 #endif
317 
318  // ******** CREATE A SINGLE STIFFNESS MATRIX, WHICH IS REPLICATED ON ALL ELEMENTS *********
319  *outStream << "Applying stiffness matrix and right hand side ... \n\n";
320 
321  // Settings and data structures for mass and stiffness matrices
323  typedef FunctionSpaceTools fst;
324  int numCells = 1;
325 
326  // Container for nodes
327  FieldContainer<double> refQuadNodes(numCells, numNodesPerElem, spaceDim);
328  // Containers for Jacobian
329  FieldContainer<double> refQuadJacobian(numCells, numCubPoints, spaceDim, spaceDim);
330  FieldContainer<double> refQuadJacobInv(numCells, numCubPoints, spaceDim, spaceDim);
331  FieldContainer<double> refQuadJacobDet(numCells, numCubPoints);
332  // Containers for element HGRAD stiffness matrix
333  FieldContainer<double> localStiffMatrix(numCells, numFieldsG, numFieldsG);
334  FieldContainer<double> weightedMeasure(numCells, numCubPoints);
335  FieldContainer<double> quadGradsTransformed(numCells, numFieldsG, numCubPoints, spaceDim);
336  FieldContainer<double> quadGradsTransformedWeighted(numCells, numFieldsG, numCubPoints, spaceDim);
337  // Containers for right hand side vectors
338  FieldContainer<double> rhsData(numCells, numCubPoints);
339  FieldContainer<double> localRHS(numCells, numFieldsG);
340  FieldContainer<double> quadGValsTransformed(numCells, numFieldsG, numCubPoints);
341  FieldContainer<double> quadGValsTransformedWeighted(numCells, numFieldsG, numCubPoints);
342  // Container for cubature points in physical space
343  FieldContainer<double> physCubPoints(numCells, numCubPoints, cubDim);
344 
345  // Global arrays in Epetra format
346  Epetra_SerialComm Comm;
347  Epetra_Map globalMapG(numDOF, 0, Comm);
348  Epetra_FEVector u(globalMapG);
349  Epetra_FEVector Ku(globalMapG);
350  u.Random();
351 
352  std::cout << "About to start ref element matrix\n";
353 
354  // ************************** Compute element HGrad stiffness matrices *******************************
355  refQuadNodes(0,0,0) = 0.0;
356  refQuadNodes(0,0,1) = 0.0;
357  refQuadNodes(0,1,0) = hx;
358  refQuadNodes(0,1,1) = 0.0;
359  refQuadNodes(0,2,0) = hx;
360  refQuadNodes(0,2,1) = hy;
361  refQuadNodes(0,3,0) = 0.0;
362  refQuadNodes(0,3,1) = hy;
363 
364  // Compute cell Jacobians, their inverses and their determinants
365  CellTools::setJacobian(refQuadJacobian, cubPoints, refQuadNodes, quad_4);
366  CellTools::setJacobianInv(refQuadJacobInv, refQuadJacobian );
367  CellTools::setJacobianDet(refQuadJacobDet, refQuadJacobian );
368 
369  // transform from [-1,1]^2 to [0,hx]x[0,hy]
370  fst::HGRADtransformGRAD<double>(quadGradsTransformed, refQuadJacobInv, quadGrads);
371 
372  // compute weighted measure
373  fst::computeCellMeasure<double>(weightedMeasure, refQuadJacobDet, cubWeights);
374 
375  // multiply values with weighted measure
376  fst::multiplyMeasure<double>(quadGradsTransformedWeighted,
377  weightedMeasure, quadGradsTransformed);
378 
379  // integrate to compute element stiffness matrix
380  fst::integrate<double>(localStiffMatrix,
381  quadGradsTransformed, quadGradsTransformedWeighted, COMP_BLAS);
382 
383  std::cout << "Finished with reference element matrix\n";
384 
385 
386  // now we will scatter global degrees of freedom, apply the local stiffness matrix
387  // with BLAS, and then gather the results
388  FieldContainer<double> uScattered(numElems,numFieldsG);
389  FieldContainer<double> KuScattered(numElems,numFieldsG);
390 
391  // to extract info from u
392 
393  u.GlobalAssemble();
394 
395  Epetra_Time multTimer(Comm);
396 
397  Ku.PutScalar(0.0);
398  Ku.GlobalAssemble();
399 
400  double *uVals = u[0];
401  double *KuVals = Ku[0];
402 
403  Teuchos::BLAS<int,double> blas;
404  Epetra_Time scatterTime(Comm);
405  std::cout << "Scattering\n";
406  // Scatter
407  for (int k=0; k<numElems; k++)
408  {
409  for (int i=0;i<numFieldsG;i++)
410  {
411  uScattered(k,i) = uVals[ltgMapping(k,i)];
412  }
413  }
414  const double scatTime = scatterTime.ElapsedTime();
415  std::cout << "Scattered in time " << scatTime << "\n";
416 
417  Epetra_Time blasTimer(Comm);
418  blas.GEMM(Teuchos::NO_TRANS , Teuchos::NO_TRANS ,
419  numFieldsG , numElems, numFieldsG ,
420  1.0 ,
421  &localStiffMatrix(0,0,0) ,
422  numFieldsG ,
423  &uScattered(0,0) ,
424  numFieldsG ,
425  0.0 ,
426  &KuScattered(0,0) ,
427  numFieldsG );
428  const double blasTime = blasTimer.ElapsedTime();
429  std::cout << "Element matrices applied in " << blasTime << "\n";
430 
431  Epetra_Time gatherTimer(Comm);
432  // Gather
433  for (int k=0;k<numElems;k++)
434  {
435  for (int i=0;i<numFieldsG;i++)
436  {
437  KuVals[ltgMapping(k,i)] += KuScattered(k,i);
438  }
439  }
440 
441  const double gatherTime = gatherTimer.ElapsedTime();
442  std::cout << "Gathered in " << gatherTime << "\n";
443 
444 
445  const double applyTime = gatherTime + blasTime + scatTime;
446  std::cout << "Time to do matrix-free product: " << applyTime << std::endl;
447 
448 
449  std::cout << "End Result: TEST PASSED\n";
450 
451  // reset format state of std::cout
452  std::cout.copyfmt(oldFormatState);
453 
454  return 0;
455 }
456 
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: