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"
68 using namespace Intrepid;
71 int main(
int argc,
char *argv[]) {
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>;
84 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
88 int iprint = argc - 1;
89 Teuchos::RCP<std::ostream> outStream;
90 Teuchos::oblackholestream bhs;
92 outStream = Teuchos::rcp(&std::cout,
false);
94 outStream = Teuchos::rcp(&bhs,
false);
96 Teuchos::oblackholestream oldFormatState;
97 oldFormatState.copyfmt(std::cout);
100 <<
"=============================================================================\n" \
102 <<
"| Unit Test (Basis_HGRAD_LINE_Hermite_FEM) |\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" \
107 <<
"| D^2[E(x) D^2 w(x) = q(x) |\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" \
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" \
117 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
118 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
119 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
121 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
122 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
124 <<
"=============================================================================\n";
131 shards::CellTopology line(shards::getCellTopologyData<shards::Line<2>>());
139 int numFields = 2*numPts;
143 double cellLength = length/numCells;
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;
155 int numDof = numCells*(numFields-2);
157 *outStream <<
"Total degrees of freedom = " << numDof << std::endl;
159 FC cellVert(numCells,numVert,spaceDim);
162 for(
int cell=0; cell<numCells; ++cell ) {
163 cellVert(cell,0,0) = cell*cellLength;
164 cellVert(cell,1,0) = (cell+1)*cellLength;
172 RCP<Cubature<double>> cellCub = cubFactory.
create(line,cubOrder);
174 FC cubPts(cubOrder, spaceDim);
176 cellCub->getCubature(cubPts,cubWts);
179 int numCubPts = cellCub->getNumPoints();
181 *outStream <<
"Number of cubature points = " << numCubPts << std::endl;
183 cubPts.resize(numCubPts,spaceDim);
184 cubWts.resize(numCubPts);
186 FC physCubPts(numCells,numCubPts, spaceDim);
187 FC wtdMeasure(numCells,numCubPts);
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);
196 *outStream << std::setw(5) << cellVert(numCells-1,1,0) << std::endl;
198 *outStream <<
"\nReference cell cubature points:" << std::endl;
199 for(
int pt=0; pt<numCubPts; ++pt) {
200 *outStream << std::setw(10) << cubPts(pt,0);
202 *outStream << std::endl;
204 *outStream <<
"\nReference cell cubature weights:" << std::endl;
205 for(
int pt=0; pt<numCubPts; ++pt) {
206 *outStream << std::setw(10) << cubWts(pt);
208 *outStream << std::endl;
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;
215 *outStream << std::endl;
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);
223 *outStream << std::endl;
231 FC elasmod(numCells,numCubPts);
232 FC qforce(numCells,numCubPts);
234 for(
int cell=0; cell<numCells; ++cell) {
236 for(
int pt=0; pt<numCubPts; ++pt) {
237 double x = physCubPts(cell,pt,0);
238 elasmod(cell,pt) = 2.0-x;
239 qforce(cell,pt) = 24.0;
248 FC pts(PointTools::getLatticeSize(line,numPts-1),1);
249 PointTools::getLattice<double,FC>(pts,line,numPts-1);
251 *outStream <<
"\nReference cell interpolation points:" << std::endl;
252 for(
int pt=0; pt<numPts; ++pt) {
253 *outStream << std::setw(10) << pts(pt,0);
255 *outStream << std::endl;
258 FC physPts(numCells, numPts, spaceDim);
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;
266 *outStream << std::endl;
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);
274 *outStream << std::endl;
279 FC valsCubPts(numFields,numCubPts);
280 FC der2CubPts(numFields,numCubPts,spaceDim);
282 hermiteBasis.getValues(valsCubPts,cubPts,OPERATOR_VALUE);
283 hermiteBasis.getValues(der2CubPts,cubPts,OPERATOR_D2);
285 FC jacobian(numCells,numCubPts,spaceDim,spaceDim);
286 FC jacInv(numCells,numCubPts,spaceDim,spaceDim);
287 FC jacDet(numCells,numCubPts);
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);
294 CT::setJacobian(jacobian,cubPts,cellVert,line);
295 CT::setJacobianInv(jacInv,jacobian);
296 CT::setJacobianDet(jacDet,jacobian);
298 FST::computeCellMeasure<double>(wtdMeasure,jacDet,cubWts);
299 FST::HGRADtransformVALUE<double>(tranValsCubPts,valsCubPts);
307 AT::matvecProductDataField<double>(tranDer2CubPts,jacInv,der2CubPts);
308 FC temp_Der2CubPts(tranDer2CubPts);
311 AT::matvecProductDataField<double>(tranDer2CubPts,jacInv,temp_Der2CubPts);
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;
330 FST::multiplyMeasure<double>(wtdTranValsCubPts,wtdMeasure,tranValsCubPts);
332 FST::multiplyMeasure<double>(wtdTranDer2CubPts,wtdMeasure,tranDer2CubPts);
333 FC temp_wtdTranDer2CubPts(wtdTranDer2CubPts);
334 FST::multiplyMeasure<double>(wtdTranDer2CubPts,elasmod,temp_wtdTranDer2CubPts);
336 FC loadVectors(numCells,numFields);
337 FC stiffnessMatrices(numCells,numFields,numFields);
339 FST::integrate(loadVectors, qforce, wtdTranValsCubPts, COMP_CPP);
340 FST::integrate(stiffnessMatrices, tranDer2CubPts, wtdTranDer2CubPts, COMP_CPP);
348 Matrix K(numDof,numDof);
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);
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);
365 for(
int cf=0; cf<numFields; ++cf ) {
366 int col = cf + (numFields-2)*cell-2;
367 K(row,col) += stiffnessMatrices(cell,rf,cf);
377 solver.setMatrix(rcpFromRef(K));
378 solver.factorWithEquilibration(
true);
380 solver.setVectors(rcpFromRef(w),rcpFromRef(q));
383 int dim = 1+numDof/2;
388 for(
int i=1; i<dim; ++i) {
394 Vector w0_exact( dim );
395 Vector w1_exact( dim );
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);
407 w0_exact(dim-1) = 1.0;
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);
417 error0 = std::sqrt(error0);
418 error1 = std::sqrt(error1);
420 *outStream <<
"\n\n";
421 *outStream <<
"|w-w_exact| = " << error0 << std::endl;
422 *outStream <<
"|w'-w'_exact| = " << error1 << std::endl;
424 double tolerance = 2e-10;
426 if( error0 > tolerance ) {
427 *outStream <<
"Solution failed to converge within tolerance" << std::endl;
431 if( error1 > tolerance ) {
432 *outStream <<
"Derivative of solution failed to converge within tolerance" << std::endl;
439 catch (
const std::logic_error & err) {
440 *outStream << err.what() <<
"\n\n";
445 std::cout <<
"End Result: TEST FAILED\n";
447 std::cout <<
"End Result: TEST PASSED\n";
450 std::cout.copyfmt(oldFormatState);
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 utility class to provide multidimensional containers.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
A factory class that generates specific instances of cubatures.
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > °ree)
Factory method.