Intrepid
test_02.cpp
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 
49 // Intrepid Includes
53 #include "Intrepid_ArrayTools.hpp"
54 #include "Intrepid_CellTools.hpp"
55 #include "Intrepid_PointTools.hpp"
57 
58 // Teuchos Includes
59 #include "Teuchos_oblackholestream.hpp"
60 #include "Teuchos_RCP.hpp"
61 #include "Teuchos_GlobalMPISession.hpp"
62 #include "Teuchos_SerialDenseSolver.hpp"
63 #include "Teuchos_SerialDenseVector.hpp"
64 
65 #include <iomanip>
66 
67 using namespace std;
68 using namespace Intrepid;
69 
70 
71 int main(int argc, char *argv[]) {
72 
73  using FC = FieldContainer<double>;
74  using CT = CellTools<double>;
75  using FST = FunctionSpaceTools;
76  using AT = ArrayTools;
77 
78  using Teuchos::RCP;
79  using Teuchos::rcpFromRef;
80  using Vector = Teuchos::SerialDenseVector<int,double>;
81  using Matrix = Teuchos::SerialDenseMatrix<int,double>;
82  using Solver = Teuchos::SerialDenseSolver<int,double>;
83 
84  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
85 
86  // This little trick lets us print to std::cout only if
87  // a (dummy) command-line argument is provided.
88  int iprint = argc - 1;
89  Teuchos::RCP<std::ostream> outStream;
90  Teuchos::oblackholestream bhs; // outputs nothing
91  if (iprint > 0)
92  outStream = Teuchos::rcp(&std::cout, false);
93  else
94  outStream = Teuchos::rcp(&bhs, false);
95  // Save the format state of the original std::cout.
96  Teuchos::oblackholestream oldFormatState;
97  oldFormatState.copyfmt(std::cout);
98 
99  *outStream \
100  << "=============================================================================\n" \
101  << "| |\n" \
102  << "| Unit Test (Basis_HGRAD_LINE_Hermite_FEM) |\n" \
103  << "| |\n" \
104  << "| Solve the cantilevered Euler-Bernoulli static beam equation with unit |\n" \
105  << "| second moment of area and verify using a manufactured solution. |\n" \
106  << "| |\n" \
107  << "| D^2[E(x) D^2 w(x) = q(x) |\n" \
108  << "| |\n" \
109  << "| with clamped boundary conditions w(0) = 0, w'(0) = 0 |\n" \
110  << "| stress boundary condition E(1)w\"(1)=-6 |\n" \
111  << "| and shear force boundary condition [Ew\"]'(1)=-6 |\n" \
112  << "| |\n" \
113  << "| The exact deflection is w(x) = 3x^2-2*x^3 |\n" \
114  << "| The elastic modulus is E(x) = 2-x |\n" \
115  << "| The forcing term is q(x) = 24 |\n" \
116  << "| |\n" \
117  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
118  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
119  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
120  << "| |\n" \
121  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
122  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
123  << "| |\n" \
124  << "=============================================================================\n";
125 
126  int errorFlag = 0;
127 
128 
129  try {
130 
131  shards::CellTopology line(shards::getCellTopologyData<shards::Line<2>>());
132 
134 
135  int numCells = 10; // Number of cells
136  int numVert = 2; // Number of vertices per cell
137  int cubOrder = 8; // Highest order of polynomial to integrate exactly
138  int numPts = 3; // Number of interpolation points per cell
139  int numFields = 2*numPts; // Number of basis functions per cell
140  int spaceDim = 1; // Number of spatial dimensions
141 
142  double length = 1.0; // Computatonal domain length
143  double cellLength = length/numCells;
144 
145  *outStream << "\n\nDiscretization Details" << std::endl;
146  *outStream << "-------------------------------" << std::endl;
147  *outStream << "Number of cells = " << numCells << std::endl;
148  *outStream << "Cubature order = " << cubOrder << std::endl;
149  *outStream << "Number of basis functions = " << numFields << std::endl;
150  *outStream << "Physical cell length = " << cellLength << std::endl;
151 
152  // Total degrees of freedom
153  // Exclude 2 for the clamped boundary condition at x=0
154  // Exclude 2 per cell for value and derivative node condensation
155  int numDof = numCells*(numFields-2);
156 
157  *outStream << "Total degrees of freedom = " << numDof << std::endl;
158 
159  FC cellVert(numCells,numVert,spaceDim); // Elemental end points
160 
161  // Set cell vertices
162  for(int cell=0; cell<numCells; ++cell ) {
163  cellVert(cell,0,0) = cell*cellLength;
164  cellVert(cell,1,0) = (cell+1)*cellLength;
165  }
166 
167 
168  /*****************************************/
169  /* CREATE CELL AND PHYSICAL CUBATURE */
170  /*****************************************/
171 
172  RCP<Cubature<double>> cellCub = cubFactory.create(line,cubOrder);
173 
174  FC cubPts(cubOrder, spaceDim); // Reference cell cubature points
175  FC cubWts(cubOrder); // Reference cell cubature weights
176  cellCub->getCubature(cubPts,cubWts);
177 
178  // Determine how many points are used and resize accordingly
179  int numCubPts = cellCub->getNumPoints();
180 
181  *outStream << "Number of cubature points = " << numCubPts << std::endl;
182 
183  cubPts.resize(numCubPts,spaceDim);
184  cubWts.resize(numCubPts);
185 
186  FC physCubPts(numCells,numCubPts, spaceDim); // Physical cubature points
187  FC wtdMeasure(numCells,numCubPts);
188 
189  CellTools<double>::mapToPhysicalFrame(physCubPts,cubPts,cellVert,line);
190 
191  *outStream << std::setprecision(5) << std::endl;
192  *outStream << "Cell Vertices:" << std::endl;;
193  for(int cell=0; cell<numCells; ++cell) {
194  *outStream << std::setw(5) << cellVert(cell,0,0);
195  }
196  *outStream << std::setw(5) << cellVert(numCells-1,1,0) << std::endl;
197 
198  *outStream << "\nReference cell cubature points:" << std::endl;
199  for(int pt=0; pt<numCubPts; ++pt) {
200  *outStream << std::setw(10) << cubPts(pt,0);
201  }
202  *outStream << std::endl;
203 
204  *outStream << "\nReference cell cubature weights:" << std::endl;
205  for( int pt=0; pt<numCubPts; ++pt) {
206  *outStream << std::setw(10) << cubWts(pt);
207  }
208  *outStream << std::endl;
209 
210  *outStream << "\nPhysical cubature points:\n" << std::endl;
211  *outStream << std::setw(7) << "Cell\\Pt | ";
212  for(int pt=0; pt<numCubPts; ++pt) {
213  *outStream << std::setw(10) << pt;
214  }
215  *outStream << std::endl;
216 
217  *outStream << std::string(10*(1+numCubPts),'-') << std::endl;
218  for(int cell=0; cell<numCells; ++cell){
219  *outStream << std::setw(7) << cell << " | ";
220  for(int pt=0; pt<numCubPts; ++pt) {
221  *outStream << std::setw(10) << physCubPts(cell,pt,0);
222  }
223  *outStream << std::endl;
224  }
225 
226 
227  /********************************************/
228  /* ELASTIC MODULUS AND FORCING FUNCTION */
229  /********************************************/
230 
231  FC elasmod(numCells,numCubPts);
232  FC qforce(numCells,numCubPts);
233 
234  for(int cell=0; cell<numCells; ++cell) {
235 
236  for(int pt=0; pt<numCubPts; ++pt) {
237  double x = physCubPts(cell,pt,0);
238  elasmod(cell,pt) = 2.0-x; //std::exp(-x);
239  qforce(cell,pt) = 24.0; // 4.0-3.0*x; //6*x; // (x-2.0)*std::exp(-x);
240  }
241  }
242 
243  /****************************************/
244  /* CREATE HERMITE INTERPOLANT BASIS */
245  /****************************************/
246 
247 
248  FC pts(PointTools::getLatticeSize(line,numPts-1),1);
249  PointTools::getLattice<double,FC>(pts,line,numPts-1);
250 
251  *outStream << "\nReference cell interpolation points:" << std::endl;
252  for(int pt=0; pt<numPts; ++pt) {
253  *outStream << std::setw(10) << pts(pt,0);
254  }
255  *outStream << std::endl;
256 
257 
258  FC physPts(numCells, numPts, spaceDim); // Physical interpolation points
259  CellTools<double>::mapToPhysicalFrame(physPts,pts,cellVert,line);
260 
261  *outStream << "\nPhysical interpolation points:\n" << std::endl;
262  *outStream << std::setw(7) << "Cell\\Pt | ";
263  for(int pt=0; pt<numPts; ++pt) {
264  *outStream << std::setw(10) << pt;
265  }
266  *outStream << std::endl;
267 
268  *outStream << std::string(10*(1+numPts),'-') << std::endl;
269  for(int cell=0; cell<numCells; ++cell){
270  *outStream << std::setw(7) << cell << " | ";
271  for(int pt=0; pt<numPts; ++pt) {
272  *outStream << std::setw(10) << physPts(cell,pt,0);
273  }
274  *outStream << std::endl;
275  }
276 
277  Basis_HGRAD_LINE_Hermite_FEM<double,FC> hermiteBasis(pts);
278 
279  FC valsCubPts(numFields,numCubPts);
280  FC der2CubPts(numFields,numCubPts,spaceDim);
281 
282  hermiteBasis.getValues(valsCubPts,cubPts,OPERATOR_VALUE);
283  hermiteBasis.getValues(der2CubPts,cubPts,OPERATOR_D2);
284 
285  FC jacobian(numCells,numCubPts,spaceDim,spaceDim);
286  FC jacInv(numCells,numCubPts,spaceDim,spaceDim);
287  FC jacDet(numCells,numCubPts);
288 
289  FC tranValsCubPts(numCells,numFields,numCubPts);
290  FC tranDer2CubPts(numCells,numFields,numCubPts,spaceDim);
291  FC wtdTranValsCubPts(numCells,numFields,numCubPts);
292  FC wtdTranDer2CubPts(numCells,numFields,numCubPts,spaceDim);
293 
294  CT::setJacobian(jacobian,cubPts,cellVert,line);
295  CT::setJacobianInv(jacInv,jacobian);
296  CT::setJacobianDet(jacDet,jacobian);
297 
298  FST::computeCellMeasure<double>(wtdMeasure,jacDet,cubWts);
299  FST::HGRADtransformVALUE<double>(tranValsCubPts,valsCubPts);
300 
301  // There is no predefined transform for second derivatives
302  // Note that applying the Jacobian inverse twice is only valid because of the
303  // affine mapping between reference and physical cells. For general nonlinear
304  // mappings, second order terms would be needed.
305 
306  // Apply once
307  AT::matvecProductDataField<double>(tranDer2CubPts,jacInv,der2CubPts);
308  FC temp_Der2CubPts(tranDer2CubPts);
309 
310  // Apply twice
311  AT::matvecProductDataField<double>(tranDer2CubPts,jacInv,temp_Der2CubPts);
312 
313 
314  // Scale derivative interpolants by cell length
315 
316  for( int cell=0; cell<numCells; ++cell ) {
317  double scale = (cellVert(cell,1,0)-cellVert(cell,0,0))/2.0;
318  for( int field=0; field<numFields/2; ++field ) {
319  for( int pt=0; pt<numCubPts; ++pt ) {
320  tranValsCubPts(cell,2*field+1,pt) *= scale;
321  tranDer2CubPts(cell,2*field+1,pt,0) *= scale;
322  }
323  }
324  }
325 
326  /********************************************/
327  /* EVALUATE FORCING AND STIFFNESS TERMS */
328  /********************************************/
329 
330  FST::multiplyMeasure<double>(wtdTranValsCubPts,wtdMeasure,tranValsCubPts);
331 
332  FST::multiplyMeasure<double>(wtdTranDer2CubPts,wtdMeasure,tranDer2CubPts);
333  FC temp_wtdTranDer2CubPts(wtdTranDer2CubPts);
334  FST::multiplyMeasure<double>(wtdTranDer2CubPts,elasmod,temp_wtdTranDer2CubPts);
335 
336  FC loadVectors(numCells,numFields);
337  FC stiffnessMatrices(numCells,numFields,numFields);
338 
339  FST::integrate(loadVectors, qforce, wtdTranValsCubPts, COMP_CPP);
340  FST::integrate(stiffnessMatrices, tranDer2CubPts, wtdTranDer2CubPts, COMP_CPP);
341 
342  /***********************************************************/
343  /* ASSEMBLY OF GLOBAL STIFFNESS MATRIX AND LOAD VECTOR */
344  /***********************************************************/
345 
346  Vector q(numDof); // Global Load Vector
347  Vector w(numDof); // Global Displacement Vector (solution)
348  Matrix K(numDof,numDof); // Global Stiffness Matrix
349 
350  // For the first cell, we exclude the first two fields to enforce the clamped
351  // boundary condition at x=0
352 
353  for( int row=0; row<numFields-2; ++row ) {
354  q(row) = loadVectors(0,row+2);
355  for(int col=0; col<numFields-2; ++col ) {
356  K(row,col) = stiffnessMatrices(0,row+2,col+2);
357  }
358  }
359 
360  for( int cell=1; cell<numCells; ++cell ) {
361  for( int rf=0; rf<numFields; ++rf ) {
362  int row = rf + (numFields-2)*cell-2;
363  q(row) += loadVectors(cell,rf);
364 
365  for( int cf=0; cf<numFields; ++cf ) {
366  int col = cf + (numFields-2)*cell-2;
367  K(row,col) += stiffnessMatrices(cell,rf,cf);
368  }
369  }
370  }
371 
372  // Boundary conditions
373  q(numDof-2) += 6.0; // Stress boundary condition
374  q(numDof-1) -= 6.0; // Shear force boundary condition
375 
376  Solver solver;
377  solver.setMatrix(rcpFromRef(K));
378  solver.factorWithEquilibration(true);
379  solver.factor();
380  solver.setVectors(rcpFromRef(w),rcpFromRef(q));
381  solver.solve();
382 
383  int dim = 1+numDof/2;
384  Vector w0( dim );
385  Vector w1( dim );
386 
387  // Separate solution into value and derivative
388  for(int i=1; i<dim; ++i) {
389  w0(i) = w(2*i-2); // Extract deflection values
390  w1(i) = w(2*i-1); // Extract deflection derivatives
391  }
392 
393  // Create exact solution and its derivative
394  Vector w0_exact( dim );
395  Vector w1_exact( dim );
396 
397  int row=0;
398  for( int cell=0; cell<numCells; ++cell ) {
399  for( int pt=0; pt<numPts-1; ++pt ) {
400  double x = physPts(cell,pt,0);
401  w0_exact(row) = (3.0-2*x)*x*x;
402  w1_exact(row) = 6.0*x*(1.0-x);
403  row++;
404  }
405  }
406 
407  w0_exact(dim-1) = 1.0;
408 
409  double error0 = 0;
410  double error1 = 0;
411 
412  for( int i=0; i<dim; ++i ) {
413  error0 += std::pow(w0(i)-w0_exact(i),2);
414  error1 += std::pow(w1(i)-w1_exact(i),2);
415  }
416 
417  error0 = std::sqrt(error0);
418  error1 = std::sqrt(error1);
419 
420  *outStream << "\n\n";
421  *outStream << "|w-w_exact| = " << error0 << std::endl;
422  *outStream << "|w'-w'_exact| = " << error1 << std::endl;
423 
424  double tolerance = 2e-10;
425 
426  if( error0 > tolerance ) {
427  *outStream << "Solution failed to converge within tolerance" << std::endl;
428  errorFlag++;
429  }
430 
431  if( error1 > tolerance ) {
432  *outStream << "Derivative of solution failed to converge within tolerance" << std::endl;
433  errorFlag++;
434  }
435 
436  }
437 
438  // Catch unexpected errors
439  catch (const std::logic_error & err) {
440  *outStream << err.what() << "\n\n";
441  errorFlag = -1000;
442  };
443 
444  if (errorFlag != 0)
445  std::cout << "End Result: TEST FAILED\n";
446  else
447  std::cout << "End Result: TEST PASSED\n";
448 
449  // reset format state of std::cout
450  std::cout.copyfmt(oldFormatState);
451 
452  return errorFlag;
453 }
Implements Hermite interpolant basis of degree n on the reference Line cell. The basis has cardinalit...
Header file for the Intrepid::HGRAD_LINE_Hermite_FEM class.
Header file for the Intrepid::CellTools class.
Header file for utility class to provide point tools, such as barycentric coordinates, equispaced lattices, and warp-blend point distrubtions.
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...
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for the Intrepid::FunctionSpaceTools class.
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays...
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: