Intrepid
test_02.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
49// Intrepid Includes
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
67using namespace std;
68using namespace Intrepid;
69
70
71int main(int argc, char *argv[]) {
72
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
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}
Header file for utility class to provide array tools, such as tensor contractions,...
Header file for the Intrepid::CellTools class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for the Intrepid::HGRAD_LINE_Hermite_FEM class.
Header file for utility class to provide point tools, such as barycentric coordinates,...
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays,...
Implements Hermite interpolant basis of degree n on the reference Line cell. The basis has cardinalit...
A stateless class for operations on cell data. Provides methods for:
static void mapToPhysicalFrame(ArrayPhysPoint &physPoints, const ArrayRefPoint &refPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computes F, the reference-to-physical frame map.
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.
Defines expert-level interfaces for the evaluation of functions and operators in physical space (supp...
static int getLatticeSize(const shards::CellTopology &cellType, const int order, const int offset=0)
Computes the number of points in a lattice of a given order on a simplex (currently disabled for othe...