Intrepid
example_03NL.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 
74 // Intrepid includes
77 #include "Intrepid_CellTools.hpp"
78 #include "Intrepid_ArrayTools.hpp"
79 #include "Intrepid_HGRAD_HEX_C1_FEM.hpp"
82 #include "Intrepid_Utils.hpp"
83 
84 // Epetra includes
85 #include "Epetra_Time.h"
86 #include "Epetra_Map.h"
87 #include "Epetra_FECrsMatrix.h"
88 #include "Epetra_FEVector.h"
89 #include "Epetra_SerialComm.h"
90 
91 // Teuchos includes
92 #include "Teuchos_oblackholestream.hpp"
93 #include "Teuchos_RCP.hpp"
94 #include "Teuchos_BLAS.hpp"
95 #include "Teuchos_Time.hpp"
96 
97 // Shards includes
98 #include "Shards_CellTopology.hpp"
99 
100 // EpetraExt includes
101 #include "EpetraExt_RowMatrixOut.h"
102 #include "EpetraExt_MultiVectorOut.h"
103 #include "EpetraExt_MatrixMatrix.h"
104 
105 // Sacado includes
106 #include "Sacado.hpp"
107 #include "Sacado_Fad_DVFad.hpp"
108 #include "Sacado_Fad_SimpleFad.hpp"
109 #include "Sacado_CacheFad_DFad.hpp"
110 #include "Sacado_CacheFad_SFad.hpp"
111 #include "Sacado_CacheFad_SLFad.hpp"
112 
113 
114 using namespace std;
115 using namespace Intrepid;
116 
117 #define INTREPID_INTEGRATE_COMP_ENGINE COMP_BLAS
118 
119 #define BATCH_SIZE 10
120 
121 //typedef Sacado::Fad::DFad<double> FadType;
122 //typedef Sacado::CacheFad::DFad<double> FadType;
123 //typedef Sacado::ELRCacheFad::DFad<double> FadType;
124 //typedef Sacado::Fad::SFad<double,8> FadType;
125 typedef Sacado::CacheFad::SFad<double,8> FadType;
126 //typedef Sacado::ELRCacheFad::SFad<double,8> FadType;
127 //typedef Sacado::Fad::SLFad<double,8> FadType;
128 //typedef Sacado::CacheFad::SLFad<double,8> FadType;
129 //typedef Sacado::ELRCacheFad::SLFad<double,8> FadType;
130 
131 //#define DUMP_DATA
132 
133 // Functions to evaluate nonlinear terms
135 
136 template<class ScalarT>
138 //
139 
140 int main(int argc, char *argv[]) {
141 
142  // Check number of arguments
143  if (argc < 4) {
144  std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
145  std::cout <<"Usage:\n\n";
146  std::cout <<" ./Intrepid_example_Drivers_Example_03NL.exe NX NY NZ verbose\n\n";
147  std::cout <<" where \n";
148  std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n";
149  std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n";
150  std::cout <<" int NZ - num intervals in z direction (assumed box domain, 0,1) \n";
151  std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
152  exit(1);
153  }
154 
155  // This little trick lets us print to std::cout only if
156  // a (dummy) command-line argument is provided.
157  int iprint = argc - 1;
158  Teuchos::RCP<std::ostream> outStream;
159  Teuchos::oblackholestream bhs; // outputs nothing
160  if (iprint > 3)
161  outStream = Teuchos::rcp(&std::cout, false);
162  else
163  outStream = Teuchos::rcp(&bhs, false);
164 
165  // Save the format state of the original std::cout.
166  Teuchos::oblackholestream oldFormatState;
167  oldFormatState.copyfmt(std::cout);
168 
169  *outStream \
170  << "===============================================================================\n" \
171  << "| |\n" \
172  << "| Example: Generate PDE Jacobian for a Nonlinear Reaction-Diffusion |\n" \
173  << "| Equation on Hexahedral Mesh |\n" \
174  << "| |\n" \
175  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
176  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
177  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
178  << "| |\n" \
179  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
180  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
181  << "| |\n" \
182  << "===============================================================================\n";
183 
184 
185  // ************************************ GET INPUTS **************************************
186 
187  int NX = atoi(argv[1]); // num intervals in x direction (assumed box domain, 0,1)
188  int NY = atoi(argv[2]); // num intervals in y direction (assumed box domain, 0,1)
189  int NZ = atoi(argv[3]); // num intervals in z direction (assumed box domain, 0,1)
190 
191  // *********************************** CELL TOPOLOGY **********************************
192 
193  // Get cell topology for base hexahedron
194  typedef shards::CellTopology CellTopology;
195  CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() );
196 
197  // Get dimensions
198  int numNodesPerElem = hex_8.getNodeCount();
199  int spaceDim = hex_8.getDimension();
200 
201  // *********************************** GENERATE MESH ************************************
202 
203  *outStream << "Generating mesh ... \n\n";
204 
205  *outStream << " NX" << " NY" << " NZ\n";
206  *outStream << std::setw(5) << NX <<
207  std::setw(5) << NY <<
208  std::setw(5) << NZ << "\n\n";
209 
210  // Print mesh information
211  int numElems = NX*NY*NZ;
212  int numNodes = (NX+1)*(NY+1)*(NZ+1);
213  *outStream << " Number of Elements: " << numElems << " \n";
214  *outStream << " Number of Nodes: " << numNodes << " \n\n";
215 
216  // Cube
217  double leftX = 0.0, rightX = 1.0;
218  double leftY = 0.0, rightY = 1.0;
219  double leftZ = 0.0, rightZ = 1.0;
220 
221  // Mesh spacing
222  double hx = (rightX-leftX)/((double)NX);
223  double hy = (rightY-leftY)/((double)NY);
224  double hz = (rightZ-leftZ)/((double)NZ);
225 
226  // Get nodal coordinates
227  FieldContainer<double> nodeCoord(numNodes, spaceDim);
228  FieldContainer<int> nodeOnBoundary(numNodes);
229  int inode = 0;
230  for (int k=0; k<NZ+1; k++) {
231  for (int j=0; j<NY+1; j++) {
232  for (int i=0; i<NX+1; i++) {
233  nodeCoord(inode,0) = leftX + (double)i*hx;
234  nodeCoord(inode,1) = leftY + (double)j*hy;
235  nodeCoord(inode,2) = leftZ + (double)k*hz;
236  if (k==0 || j==0 || i==0 || k==NZ || j==NY || i==NX){
237  nodeOnBoundary(inode)=1;
238  }
239  else {
240  nodeOnBoundary(inode)=0;
241  }
242  inode++;
243  }
244  }
245  }
246 
247 #ifdef DUMP_DATA
248  // Print nodal coords
249  ofstream fcoordout("coords.dat");
250  for (int i=0; i<numNodes; i++) {
251  fcoordout << nodeCoord(i,0) <<" ";
252  fcoordout << nodeCoord(i,1) <<" ";
253  fcoordout << nodeCoord(i,2) <<"\n";
254  }
255  fcoordout.close();
256 #endif
257 
258 
259  // Element to Node map
260  FieldContainer<int> elemToNode(numElems, numNodesPerElem);
261  int ielem = 0;
262  for (int k=0; k<NZ; k++) {
263  for (int j=0; j<NY; j++) {
264  for (int i=0; i<NX; i++) {
265  elemToNode(ielem,0) = (NY + 1)*(NX + 1)*k + (NX + 1)*j + i;
266  elemToNode(ielem,1) = (NY + 1)*(NX + 1)*k + (NX + 1)*j + i + 1;
267  elemToNode(ielem,2) = (NY + 1)*(NX + 1)*k + (NX + 1)*(j + 1) + i + 1;
268  elemToNode(ielem,3) = (NY + 1)*(NX + 1)*k + (NX + 1)*(j + 1) + i;
269  elemToNode(ielem,4) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*j + i;
270  elemToNode(ielem,5) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*j + i + 1;
271  elemToNode(ielem,6) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*(j + 1) + i + 1;
272  elemToNode(ielem,7) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*(j + 1) + i;
273  ielem++;
274  }
275  }
276  }
277 #ifdef DUMP_DATA
278  // Output connectivity
279  ofstream fe2nout("elem2node.dat");
280  for (int k=0; k<NZ; k++) {
281  for (int j=0; j<NY; j++) {
282  for (int i=0; i<NX; i++) {
283  int ielem = i + j * NX + k * NX * NY;
284  for (int m=0; m<numNodesPerElem; m++){
285  fe2nout << elemToNode(ielem,m) <<" ";
286  }
287  fe2nout <<"\n";
288  }
289  }
290  }
291  fe2nout.close();
292 #endif
293 
294 
295  // ************************************ CUBATURE **************************************
296 
297  *outStream << "Getting cubature ... \n\n";
298 
299  // Get numerical integration points and weights
300  DefaultCubatureFactory<double> cubFactory;
301  int cubDegree = 2;
302  Teuchos::RCP<Cubature<double> > hexCub = cubFactory.create(hex_8, cubDegree);
303 
304  int cubDim = hexCub->getDimension();
305  int numCubPoints = hexCub->getNumPoints();
306 
307  FieldContainer<double> cubPoints(numCubPoints, cubDim);
308  FieldContainer<double> cubWeights(numCubPoints);
309 
310  hexCub->getCubature(cubPoints, cubWeights);
311 
312 
313  // ************************************** BASIS ***************************************
314 
315  *outStream << "Getting basis ... \n\n";
316 
317  // Define basis
319  int numFieldsG = hexHGradBasis.getCardinality();
320  FieldContainer<double> hexGVals(numFieldsG, numCubPoints);
321  FieldContainer<double> hexGrads(numFieldsG, numCubPoints, spaceDim);
322 
323  // Evaluate basis values and gradients at cubature points
324  hexHGradBasis.getValues(hexGVals, cubPoints, OPERATOR_VALUE);
325  hexHGradBasis.getValues(hexGrads, cubPoints, OPERATOR_GRAD);
326 
327 
328  // ******** FEM ASSEMBLY *************
329 
330  *outStream << "Building PDE Jacobian ... \n\n";
331 
332  // Settings and data structures for mass and stiffness matrices
334  typedef FunctionSpaceTools fst;
335  int numCells = BATCH_SIZE;
336  int numBatches = numElems/numCells;
337 
338  // Container for nodes
339  FieldContainer<double> hexNodes(numCells, numNodesPerElem, spaceDim);
340  // Containers for Jacobian
341  FieldContainer<double> hexJacobian(numCells, numCubPoints, spaceDim, spaceDim);
342  FieldContainer<double> hexJacobInv(numCells, numCubPoints, spaceDim, spaceDim);
343  FieldContainer<double> hexJacobDet(numCells, numCubPoints);
344  // Containers for HGRAD bases
345  FieldContainer<double> localPDEjacobian(numCells, numFieldsG, numFieldsG);
346  FieldContainer<double> weightedMeasure(numCells, numCubPoints);
347  FieldContainer<double> hexGValsTransformed(numCells, numFieldsG, numCubPoints);
348  FieldContainer<double> hexGValsTransformedWeighted(numCells, numFieldsG, numCubPoints);
349  FieldContainer<double> hexGradsTransformed(numCells, numFieldsG, numCubPoints, spaceDim);
350  FieldContainer<double> hexGradsTransformedWeighted(numCells, numFieldsG, numCubPoints, spaceDim);
351 
352  // Global arrays in Epetra format
353  Epetra_SerialComm Comm;
354  Epetra_Map globalMapG(numNodes, 0, Comm);
355  Epetra_FECrsMatrix StiffMatrix(Copy, globalMapG, 64);
356 
357  // Additional arrays used in analytic assembly
358  FieldContainer<double> u_coeffs(numCells, numFieldsG);
359  FieldContainer<double> u_FE_val(numCells, numCubPoints);
360  FieldContainer<double> df_of_u(numCells, numCubPoints);
361  FieldContainer<double> df_of_u_times_basis(numCells, numFieldsG, numCubPoints);
362 
363 
364  // Additional arrays used in AD-based assembly.
365  FieldContainer<FadType> u_coeffsAD(numCells, numFieldsG);
366  FieldContainer<FadType> u_FE_gradAD(numCells, numCubPoints, spaceDim);
367  FieldContainer<FadType> u_FE_valAD(numCells, numCubPoints);
368  FieldContainer<FadType> f_of_u_AD(numCells, numCubPoints);
369  FieldContainer<FadType> cellResidualAD(numCells, numFieldsG);
370  for (int c=0; c<numCells; c++) {
371  for(int f=0; f<numFieldsG; f++) {
372  u_coeffsAD(c,f) = FadType(numFieldsG, f, 1.3);
373  }
374  }
375 
376  Teuchos::Time timer_jac_analytic("Time to compute element PDE Jacobians analytically: ");
377  Teuchos::Time timer_jac_fad ("Time to compute element PDE Jacobians using AD: ");
378  Teuchos::Time timer_jac_insert ("Time for global insert, w/o graph: ");
379  Teuchos::Time timer_jac_insert_g("Time for global insert, w/ graph: ");
380  Teuchos::Time timer_jac_ga ("Time for GlobalAssemble, w/o graph: ");
381  Teuchos::Time timer_jac_ga_g ("Time for GlobalAssemble, w/ graph: ");
382  Teuchos::Time timer_jac_fc ("Time for FillComplete, w/o graph: ");
383  Teuchos::Time timer_jac_fc_g ("Time for FillComplete, w/ graph: ");
384 
385 
386 
387 
388  // *** Analytic element loop ***
389  for (int bi=0; bi<numBatches; bi++) {
390 
391  // Physical cell coordinates
392  for (int ci=0; ci<numCells; ci++) {
393  int k = bi*numCells+ci;
394  for (int i=0; i<numNodesPerElem; i++) {
395  hexNodes(ci,i,0) = nodeCoord(elemToNode(k,i),0);
396  hexNodes(ci,i,1) = nodeCoord(elemToNode(k,i),1);
397  hexNodes(ci,i,2) = nodeCoord(elemToNode(k,i),2);
398  }
399  }
400 
401  // Compute cell Jacobians, their inverses and their determinants
402  CellTools::setJacobian(hexJacobian, cubPoints, hexNodes, hex_8);
403  CellTools::setJacobianInv(hexJacobInv, hexJacobian );
404  CellTools::setJacobianDet(hexJacobDet, hexJacobian );
405 
406  // ******************** COMPUTE ELEMENT HGrad STIFFNESS MATRICES WITHOUT AD *******************
407 
408  // transform to physical coordinates
409  fst::HGRADtransformGRAD<double>(hexGradsTransformed, hexJacobInv, hexGrads);
410 
411  // compute weighted measure
412  fst::computeCellMeasure<double>(weightedMeasure, hexJacobDet, cubWeights);
413 
414  // multiply values with weighted measure
415  fst::multiplyMeasure<double>(hexGradsTransformedWeighted,
416  weightedMeasure, hexGradsTransformed);
417 
418  // u_coeffs equals the value of u_coeffsAD
419  for(int i=0; i<numFieldsG; i++){
420  u_coeffs(0,i) = u_coeffsAD(0,i).val();
421  }
422 
423  timer_jac_analytic.start(); // START TIMER
424  // integrate to account for linear stiffness term
425  fst::integrate<double>(localPDEjacobian, hexGradsTransformed, hexGradsTransformedWeighted, INTREPID_INTEGRATE_COMP_ENGINE);
426 
427  // represent value of the current state (iterate) as a linear combination of the basis functions
428  u_FE_val.initialize();
429  fst::evaluate<double>(u_FE_val, u_coeffs, hexGValsTransformed);
430 
431  // evaluate derivative of the nonlinear term and multiply by basis function
432  dfunc_u(df_of_u, u_FE_val);
433  fst::scalarMultiplyDataField<double>(df_of_u_times_basis, df_of_u, hexGValsTransformed);
434 
435  // integrate to account for nonlinear reaction term
436  fst::integrate<double>(localPDEjacobian, df_of_u_times_basis, hexGValsTransformedWeighted, INTREPID_INTEGRATE_COMP_ENGINE, true);
437  timer_jac_analytic.stop(); // STOP TIMER
438 
439  // assemble into global matrix
440  for (int ci=0; ci<numCells; ci++) {
441  int k = bi*numCells+ci;
442  std::vector<int> rowIndex(numFieldsG);
443  std::vector<int> colIndex(numFieldsG);
444  for (int row = 0; row < numFieldsG; row++){
445  rowIndex[row] = elemToNode(k,row);
446  }
447  for (int col = 0; col < numFieldsG; col++){
448  colIndex[col] = elemToNode(k,col);
449  }
450  // We can insert an entire matrix at a time, but we opt for rows only.
451  //timer_jac_insert.start();
452  //StiffMatrix.InsertGlobalValues(numFieldsG, &rowIndex[0], numFieldsG, &colIndex[0], &localPDEjacobian(ci,0,0));
453  //timer_jac_insert.stop();
454  for (int row = 0; row < numFieldsG; row++){
455  timer_jac_insert.start();
456  StiffMatrix.InsertGlobalValues(1, &rowIndex[row], numFieldsG, &colIndex[0], &localPDEjacobian(ci,row,0));
457  timer_jac_insert.stop();
458  }
459  }
460 
461  } // *** end analytic element loop ***
462 
463  // Assemble global objects
464  timer_jac_ga.start(); StiffMatrix.GlobalAssemble(); timer_jac_ga.stop();
465  timer_jac_fc.start(); StiffMatrix.FillComplete(); timer_jac_fc.stop();
466 
467 
468 
469 
470  // *** AD element loop ***
471 
472  Epetra_CrsGraph mgraph = StiffMatrix.Graph();
473  Epetra_FECrsMatrix StiffMatrixViaAD(Copy, mgraph);
474 
475  for (int bi=0; bi<numBatches; bi++) {
476 
477  // ******************** COMPUTE ELEMENT HGrad STIFFNESS MATRICES AND RIGHT-HAND SIDE WITH AD ********************
478 
479  // Physical cell coordinates
480  for (int ci=0; ci<numCells; ci++) {
481  int k = bi*numCells+ci;
482  for (int i=0; i<numNodesPerElem; i++) {
483  hexNodes(ci,i,0) = nodeCoord(elemToNode(k,i),0);
484  hexNodes(ci,i,1) = nodeCoord(elemToNode(k,i),1);
485  hexNodes(ci,i,2) = nodeCoord(elemToNode(k,i),2);
486  }
487  }
488 
489  // Compute cell Jacobians, their inverses and their determinants
490  CellTools::setJacobian(hexJacobian, cubPoints, hexNodes, hex_8);
491  CellTools::setJacobianInv(hexJacobInv, hexJacobian );
492  CellTools::setJacobianDet(hexJacobDet, hexJacobian );
493 
494  // transform to physical coordinates
495  fst::HGRADtransformGRAD<double>(hexGradsTransformed, hexJacobInv, hexGrads);
496 
497  // compute weighted measure
498  fst::computeCellMeasure<double>(weightedMeasure, hexJacobDet, cubWeights);
499 
500  // multiply values with weighted measure
501  fst::multiplyMeasure<double>(hexGradsTransformedWeighted, weightedMeasure, hexGradsTransformed);
502 
503  // transform basis values to physical coordinates
504  fst::HGRADtransformVALUE<double>(hexGValsTransformed, hexGVals);
505 
506  // multiply values with weighted measure
507  fst::multiplyMeasure<double>(hexGValsTransformedWeighted,
508  weightedMeasure, hexGValsTransformed);
509 
510  timer_jac_fad.start(); // START TIMER
511  // represent gradient of the current state (iterate) as a linear combination of the gradients of basis functions
512  // use AD arrays !
513  u_FE_gradAD.initialize();
514  fst::evaluate<FadType>(u_FE_gradAD, u_coeffsAD, hexGradsTransformed);
515 
516  // represent value of the current state (iterate) as a linear combination of the basis functions
517  // use AD arrays !
518  u_FE_valAD.initialize();
519  fst::evaluate<FadType>(u_FE_valAD, u_coeffsAD, hexGValsTransformed);
520  // compute nonlinear term
521  func_u(f_of_u_AD, u_FE_valAD);
522 
523  // integrate to compute element residual
524  fst::integrate<FadType>(cellResidualAD, u_FE_gradAD, hexGradsTransformedWeighted, INTREPID_INTEGRATE_COMP_ENGINE);
525  fst::integrate<FadType>(cellResidualAD, f_of_u_AD, hexGValsTransformedWeighted, INTREPID_INTEGRATE_COMP_ENGINE, true);
526  timer_jac_fad.stop(); // STOP TIMER
527 
528  // assemble into global matrix
529  for (int ci=0; ci<numCells; ci++) {
530  int k = bi*numCells+ci;
531  std::vector<int> rowIndex(numFieldsG);
532  std::vector<int> colIndex(numFieldsG);
533  for (int row = 0; row < numFieldsG; row++){
534  rowIndex[row] = elemToNode(k,row);
535  }
536  for (int col = 0; col < numFieldsG; col++){
537  colIndex[col] = elemToNode(k,col);
538  }
539  for (int row = 0; row < numFieldsG; row++){
540  timer_jac_insert_g.start();
541  StiffMatrixViaAD.SumIntoGlobalValues(1, &rowIndex[row], numFieldsG, &colIndex[0], cellResidualAD(ci,row).dx());
542  timer_jac_insert_g.stop();
543  }
544  }
545 
546  } // *** end AD element loop ***
547 
548  // Assemble global objects
549  timer_jac_ga_g.start(); StiffMatrixViaAD.GlobalAssemble(); timer_jac_ga_g.stop();
550  timer_jac_fc_g.start(); StiffMatrixViaAD.FillComplete(); timer_jac_fc_g.stop();
551 
552 
553 
554  /****** Output *******/
555 
556 #ifdef DUMP_DATA
557  // Dump matrices to disk
558  EpetraExt::RowMatrixToMatlabFile("stiff_matrix.dat",StiffMatrix);
559  EpetraExt::RowMatrixToMatlabFile("stiff_matrixAD.dat",StiffMatrixViaAD);
560 #endif
561 
562  // take the infinity norm of the difference between StiffMatrix and StiffMatrixViaAD to see that
563  // the two matrices are the same
564  EpetraExt::MatrixMatrix::Add(StiffMatrix, false, 1.0, StiffMatrixViaAD, -1.0);
565  double normMat = StiffMatrixViaAD.NormInf();
566  *outStream << "Infinity norm of difference between stiffness matrices = " << normMat << "\n";
567 
568 
569  *outStream << "\n\nNumber of global nonzeros: " << StiffMatrix.NumGlobalNonzeros() << "\n\n";
570 
571  *outStream << timer_jac_analytic.name() << " " << timer_jac_analytic.totalElapsedTime() << " sec\n";
572  *outStream << timer_jac_fad.name() << " " << timer_jac_fad.totalElapsedTime() << " sec\n\n";
573  *outStream << timer_jac_insert.name() << " " << timer_jac_insert.totalElapsedTime() << " sec\n";
574  *outStream << timer_jac_insert_g.name() << " " << timer_jac_insert_g.totalElapsedTime() << " sec\n\n";
575  *outStream << timer_jac_ga.name() << " " << timer_jac_ga.totalElapsedTime() << " sec\n";
576  *outStream << timer_jac_ga_g.name() << " " << timer_jac_ga_g.totalElapsedTime() << " sec\n\n";
577  *outStream << timer_jac_fc.name() << " " << timer_jac_fc.totalElapsedTime() << " sec\n";
578  *outStream << timer_jac_fc_g.name() << " " << timer_jac_fc_g.totalElapsedTime() << " sec\n\n";
579 
580  if ((normMat < 1.0e4*INTREPID_TOL)) {
581  std::cout << "End Result: TEST PASSED\n";
582  }
583  else {
584  std::cout << "End Result: TEST FAILED\n";
585  }
586 
587  // reset format state of std::cout
588  std::cout.copyfmt(oldFormatState);
589 
590  return 0;
591 }
592 
593 
594 template<class ScalarT>
596  int num_cells = u.dimension(0);
597  int num_cub_p = u.dimension(1);
598  for(int c=0; c<num_cells; c++){
599  for(int p=0; p<num_cub_p; p++){
600  fu(c,p) = std::pow(u(c,p),3) + std::exp(u(c,p));
601  }
602  }
603 }
604 
605 
606 void dfunc_u(FieldContainer<double> dfu, FieldContainer<double> u) {
607  int num_cells = u.dimension(0);
608  int num_cub_p = u.dimension(1);
609  for(int c=0; c<num_cells; c++) {
610  for(int p=0; p<num_cub_p; p++) {
611  dfu(c,p) = 3*u(c,p)*u(c,p) + std::exp(u(c,p));
612  }
613  }
614 }
virtual int getCardinality() const
Returns cardinality of the basis.
Header file for the Intrepid::CellTools class.
int dimension(const int whichDim) const
Returns the specified dimension.
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...
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Hexahedron cell...
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.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
A stateless class for operations on cell data. Provides methods for: