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
57#include "Teuchos_oblackholestream.hpp"
58#include "Teuchos_RCP.hpp"
59#include "Teuchos_GlobalMPISession.hpp"
60#include "Teuchos_SerialDenseMatrix.hpp"
61#include "Teuchos_SerialDenseVector.hpp"
62#include "Teuchos_LAPACK.hpp"
63
64using namespace std;
65using namespace Intrepid;
66
67void rhsFunc(FieldContainer<double> &, const FieldContainer<double> &, int, int, int);
71 const shards::CellTopology & ,
72 int, int, int, int);
73void u_exact(FieldContainer<double> &, const FieldContainer<double> &, int, int, int);
74
77 const FieldContainer<double> & points,
78 int xd,
79 int yd,
80 int zd) {
81
82 int x = 0, y = 1, z = 2;
83
84 // second x-derivatives of u
85 if (xd > 1) {
86 for (int cell=0; cell<result.dimension(0); cell++) {
87 for (int pt=0; pt<result.dimension(1); pt++) {
88 result(cell,pt) = - xd*(xd-1)*std::pow(points(cell,pt,x), xd-2) *
89 std::pow(points(cell,pt,y), yd) * std::pow(points(cell,pt,z), zd);
90 }
91 }
92 }
93
94 // second y-derivatives of u
95 if (yd > 1) {
96 for (int cell=0; cell<result.dimension(0); cell++) {
97 for (int pt=0; pt<result.dimension(1); pt++) {
98 result(cell,pt) -= yd*(yd-1)*std::pow(points(cell,pt,y), yd-2) *
99 std::pow(points(cell,pt,x), xd) * std::pow(points(cell,pt,z), zd);
100 }
101 }
102 }
103
104 // second z-derivatives of u
105 if (zd > 1) {
106 for (int cell=0; cell<result.dimension(0); cell++) {
107 for (int pt=0; pt<result.dimension(1); pt++) {
108 result(cell,pt) -= zd*(zd-1)*std::pow(points(cell,pt,z), zd-2) *
109 std::pow(points(cell,pt,x), xd) * std::pow(points(cell,pt,y), yd);
110 }
111 }
112 }
113
114 // add u
115 for (int cell=0; cell<result.dimension(0); cell++) {
116 for (int pt=0; pt<result.dimension(1); pt++) {
117 result(cell,pt) += std::pow(points(cell,pt,x), xd) * std::pow(points(cell,pt,y), yd) * std::pow(points(cell,pt,z), zd);
118 }
119 }
120
121}
122
123
126 const FieldContainer<double> & points,
127 const FieldContainer<double> & jacs,
128 const shards::CellTopology & parentCell,
129 int sideOrdinal, int xd, int yd, int zd) {
130
131 int x = 0, y = 1, z = 2;
132
133 int numCells = result.dimension(0);
134 int numPoints = result.dimension(1);
135
136 FieldContainer<double> grad_u(numCells, numPoints, 3);
137 FieldContainer<double> side_normals(numCells, numPoints, 3);
138 FieldContainer<double> normal_lengths(numCells, numPoints);
139
140 // first x-derivatives of u
141 if (xd > 0) {
142 for (int cell=0; cell<numCells; cell++) {
143 for (int pt=0; pt<numPoints; pt++) {
144 grad_u(cell,pt,x) = xd*std::pow(points(cell,pt,x), xd-1) *
145 std::pow(points(cell,pt,y), yd) * std::pow(points(cell,pt,z), zd);
146 }
147 }
148 }
149
150 // first y-derivatives of u
151 if (yd > 0) {
152 for (int cell=0; cell<numCells; cell++) {
153 for (int pt=0; pt<numPoints; pt++) {
154 grad_u(cell,pt,y) = yd*std::pow(points(cell,pt,y), yd-1) *
155 std::pow(points(cell,pt,x), xd) * std::pow(points(cell,pt,z), zd);
156 }
157 }
158 }
159
160 // first z-derivatives of u
161 if (zd > 0) {
162 for (int cell=0; cell<numCells; cell++) {
163 for (int pt=0; pt<numPoints; pt++) {
164 grad_u(cell,pt,z) = zd*std::pow(points(cell,pt,z), zd-1) *
165 std::pow(points(cell,pt,x), xd) * std::pow(points(cell,pt,y), yd);
166 }
167 }
168 }
169
170 CellTools<double>::getPhysicalSideNormals(side_normals, jacs, sideOrdinal, parentCell);
171
172 // scale normals
173 RealSpaceTools<double>::vectorNorm(normal_lengths, side_normals, NORM_TWO);
174 FunctionSpaceTools::scalarMultiplyDataData<double>(side_normals, normal_lengths, side_normals, true);
175
176 FunctionSpaceTools::dotMultiplyDataData<double>(result, grad_u, side_normals);
177
178}
179
181void u_exact(FieldContainer<double> & result, const FieldContainer<double> & points, int xd, int yd, int zd) {
182 int x = 0, y = 1, z = 2;
183 for (int cell=0; cell<result.dimension(0); cell++) {
184 for (int pt=0; pt<result.dimension(1); pt++) {
185 result(cell,pt) = std::pow(points(pt,x), xd)*std::pow(points(pt,y), yd)*std::pow(points(pt,z), zd);
186 }
187 }
188}
189
190
191
192
193int main(int argc, char *argv[]) {
194
195 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
196
197 // This little trick lets us print to std::cout only if
198 // a (dummy) command-line argument is provided.
199 int iprint = argc - 1;
200 Teuchos::RCP<std::ostream> outStream;
201 Teuchos::oblackholestream bhs; // outputs nothing
202 if (iprint > 0)
203 outStream = Teuchos::rcp(&std::cout, false);
204 else
205 outStream = Teuchos::rcp(&bhs, false);
206
207 // Save the format state of the original std::cout.
208 Teuchos::oblackholestream oldFormatState;
209 oldFormatState.copyfmt(std::cout);
210
211 *outStream \
212 << "===============================================================================\n" \
213 << "| |\n" \
214 << "| Unit Test (Basis_HGRAD_HEX_Cn_FEM) |\n" \
215 << "| |\n" \
216 << "| 1) Patch test involving mass and stiffness matrices, |\n" \
217 << "| for the Neumann problem on a physical parallelepiped |\n" \
218 << "| AND a reference hex Omega with boundary Gamma. |\n" \
219 << "| |\n" \
220 << "| - div (grad u) + u = f in Omega, (grad u) . n = g on Gamma |\n" \
221 << "| |\n" \
222 << "| For a generic parallelepiped, the basis recovers a complete |\n" \
223 << "| polynomial space of order 2. On a (scaled and/or translated) |\n" \
224 << "| reference hex, the basis recovers a complete tensor product |\n" \
225 << "| space of order 1 (i.e. incl. cross terms, e.g. x^2*y^2*z^2). |\n" \
226 << "| |\n" \
227 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
228 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
229 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
230 << "| |\n" \
231 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
232 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
233 << "| |\n" \
234 << "===============================================================================\n"\
235 << "| TEST 1: Patch test |\n"\
236 << "===============================================================================\n";
237
238
239 int errorFlag = 0;
240
241 outStream -> precision(16);
242
243
244 try {
245
246 int max_order = 3; // max total order of polynomial solution
247 DefaultCubatureFactory<double> cubFactory; // create factory
248 shards::CellTopology cell(shards::getCellTopologyData< shards::Hexahedron<> >()); // create parent cell topology
249 shards::CellTopology side(shards::getCellTopologyData< shards::Quadrilateral<> >()); // create relevant subcell (side) topology
250 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >()); // create relevant subcell (side) topology
251 int cellDim = cell.getDimension();
252 int sideDim = side.getDimension();
253 unsigned numSides = 6;
254
255 // Define array containing points at which the solution is evaluated, on the reference tet.
256 int numIntervals = 10;
257 int numInterpPoints = (numIntervals + 1)*(numIntervals + 1)*(numIntervals + 1);
258 FieldContainer<double> interp_points_ref(numInterpPoints, 3);
259 int counter = 0;
260 for (int k=0; k<=numIntervals; k++) {
261 for (int j=0; j<=numIntervals; j++) {
262 for (int i=0; i<=numIntervals; i++) {
263 interp_points_ref(counter,0) = i*(1.0/numIntervals)-1.0;
264 interp_points_ref(counter,1) = j*(1.0/numIntervals)-1.0;
265 interp_points_ref(counter,2) = k*(1.0/numIntervals)-1.0;
266 counter++;
267 }
268 }
269 }
270
271 /* Parent cell definition. */
272 FieldContainer<double> cell_nodes[2];
273 cell_nodes[0].resize(1, 8, cellDim);
274 cell_nodes[1].resize(1, 8, cellDim);
275
276 // Generic parallelepiped.
277 cell_nodes[0](0, 0, 0) = -5.0;
278 cell_nodes[0](0, 0, 1) = -1.0;
279 cell_nodes[0](0, 0, 2) = 0.0;
280 cell_nodes[0](0, 1, 0) = 4.0;
281 cell_nodes[0](0, 1, 1) = 1.0;
282 cell_nodes[0](0, 1, 2) = 1.0;
283 cell_nodes[0](0, 2, 0) = 8.0;
284 cell_nodes[0](0, 2, 1) = 3.0;
285 cell_nodes[0](0, 2, 2) = 1.0;
286 cell_nodes[0](0, 3, 0) = -1.0;
287 cell_nodes[0](0, 3, 1) = 1.0;
288 cell_nodes[0](0, 3, 2) = 0.0;
289 cell_nodes[0](0, 4, 0) = 5.0;
290 cell_nodes[0](0, 4, 1) = 9.0;
291 cell_nodes[0](0, 4, 2) = 1.0;
292 cell_nodes[0](0, 5, 0) = 14.0;
293 cell_nodes[0](0, 5, 1) = 11.0;
294 cell_nodes[0](0, 5, 2) = 2.0;
295 cell_nodes[0](0, 6, 0) = 18.0;
296 cell_nodes[0](0, 6, 1) = 13.0;
297 cell_nodes[0](0, 6, 2) = 2.0;
298 cell_nodes[0](0, 7, 0) = 9.0;
299 cell_nodes[0](0, 7, 1) = 11.0;
300 cell_nodes[0](0, 7, 2) = 1.0;
301 // Reference hex.
302 cell_nodes[1](0, 0, 0) = -1.0;
303 cell_nodes[1](0, 0, 1) = -1.0;
304 cell_nodes[1](0, 0, 2) = -1.0;
305 cell_nodes[1](0, 1, 0) = 1.0;
306 cell_nodes[1](0, 1, 1) = -1.0;
307 cell_nodes[1](0, 1, 2) = -1.0;
308 cell_nodes[1](0, 2, 0) = 1.0;
309 cell_nodes[1](0, 2, 1) = 1.0;
310 cell_nodes[1](0, 2, 2) = -1.0;
311 cell_nodes[1](0, 3, 0) = -1.0;
312 cell_nodes[1](0, 3, 1) = 1.0;
313 cell_nodes[1](0, 3, 2) = -1.0;
314 cell_nodes[1](0, 4, 0) = -1.0;
315 cell_nodes[1](0, 4, 1) = -1.0;
316 cell_nodes[1](0, 4, 2) = 1.0;
317 cell_nodes[1](0, 5, 0) = 1.0;
318 cell_nodes[1](0, 5, 1) = -1.0;
319 cell_nodes[1](0, 5, 2) = 1.0;
320 cell_nodes[1](0, 6, 0) = 1.0;
321 cell_nodes[1](0, 6, 1) = 1.0;
322 cell_nodes[1](0, 6, 2) = 1.0;
323 cell_nodes[1](0, 7, 0) = -1.0;
324 cell_nodes[1](0, 7, 1) = 1.0;
325 cell_nodes[1](0, 7, 2) = 1.0;
326
327 std::stringstream mystream[2];
328 mystream[0].str("\n>> Now testing basis on a generic parallelepiped ...\n");
329 mystream[1].str("\n>> Now testing basis on the reference hex ...\n");
330
331
332 for (int pcell = 0; pcell < 2; pcell++) {
333 *outStream << mystream[pcell].str();
334 FieldContainer<double> interp_points(1, numInterpPoints, cellDim);
335 CellTools<double>::mapToPhysicalFrame(interp_points, interp_points_ref, cell_nodes[pcell], cell);
336 interp_points.resize(numInterpPoints, cellDim);
337
338 for (int x_order=0; x_order <= max_order; x_order++) {
339 int max_y_order = max_order;
340 if (pcell == 0) {
341 max_y_order -= x_order;
342 }
343 for (int y_order=0; y_order <= max_y_order; y_order++) {
344 int max_z_order = max_order;
345 if (pcell == 0) {
346 max_z_order -= x_order;
347 max_z_order -= y_order;
348 }
349 for (int z_order=0; z_order <= max_z_order; z_order++) {
350
351 // evaluate exact solution
352 FieldContainer<double> exact_solution(1, numInterpPoints);
353 u_exact(exact_solution, interp_points, x_order, y_order, z_order);
354
355 int basis_order = max_order;
356
357 // set test tolerance;
358 double zero = basis_order*basis_order*basis_order*100*INTREPID_TOL;
359
360 //create basis
361 FieldContainer<double> pts(PointTools::getLatticeSize(line,basis_order),1);
362 PointTools::getLattice<double,FieldContainer<double> >(pts,line,basis_order);
363
364 Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
365 Teuchos::rcp(new Basis_HGRAD_HEX_Cn_FEM<double,FieldContainer<double> >( basis_order, POINTTYPE_SPECTRAL ) );
366 int numFields = basis->getCardinality();
367
368 // create cubatures
369 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*basis_order);
370 Teuchos::RCP<Cubature<double> > sideCub = cubFactory.create(side, 2*basis_order);
371 int numCubPointsCell = cellCub->getNumPoints();
372 int numCubPointsSide = sideCub->getNumPoints();
373
374 /* Computational arrays. */
375 /* Section 1: Related to parent cell integration. */
376 FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
377 FieldContainer<double> cub_points_cell_physical(1, numCubPointsCell, cellDim);
378 FieldContainer<double> cub_weights_cell(numCubPointsCell);
379 FieldContainer<double> jacobian_cell(1, numCubPointsCell, cellDim, cellDim);
380 FieldContainer<double> jacobian_inv_cell(1, numCubPointsCell, cellDim, cellDim);
381 FieldContainer<double> jacobian_det_cell(1, numCubPointsCell);
382 FieldContainer<double> weighted_measure_cell(1, numCubPointsCell);
383
384 FieldContainer<double> value_of_basis_at_cub_points_cell(numFields, numCubPointsCell);
385 FieldContainer<double> transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
386 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
387 FieldContainer<double> grad_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim);
388 FieldContainer<double> transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
389 FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
390 FieldContainer<double> fe_matrix(1, numFields, numFields);
391
392 FieldContainer<double> rhs_at_cub_points_cell_physical(1, numCubPointsCell);
393 FieldContainer<double> rhs_and_soln_vector(1, numFields);
394
395 /* Section 2: Related to subcell (side) integration. */
396 FieldContainer<double> cub_points_side(numCubPointsSide, sideDim);
397 FieldContainer<double> cub_weights_side(numCubPointsSide);
398 FieldContainer<double> cub_points_side_refcell(numCubPointsSide, cellDim);
399 FieldContainer<double> cub_points_side_physical(1, numCubPointsSide, cellDim);
400 FieldContainer<double> jacobian_side_refcell(1, numCubPointsSide, cellDim, cellDim);
401 FieldContainer<double> jacobian_det_side_refcell(1, numCubPointsSide);
402 FieldContainer<double> weighted_measure_side_refcell(1, numCubPointsSide);
403
404 FieldContainer<double> value_of_basis_at_cub_points_side_refcell(numFields, numCubPointsSide);
405 FieldContainer<double> transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
406 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
407 FieldContainer<double> neumann_data_at_cub_points_side_physical(1, numCubPointsSide);
408 FieldContainer<double> neumann_fields_per_side(1, numFields);
409
410 /* Section 3: Related to global interpolant. */
411 FieldContainer<double> value_of_basis_at_interp_points_ref(numFields, numInterpPoints);
412 FieldContainer<double> transformed_value_of_basis_at_interp_points_ref(1, numFields, numInterpPoints);
413 FieldContainer<double> interpolant(1, numInterpPoints);
414
415 FieldContainer<int> ipiv(numFields);
416
417
418
419 /******************* START COMPUTATION ***********************/
420
421 // get cubature points and weights
422 cellCub->getCubature(cub_points_cell, cub_weights_cell);
423
424 // compute geometric cell information
425 CellTools<double>::setJacobian(jacobian_cell, cub_points_cell, cell_nodes[pcell], cell);
426 CellTools<double>::setJacobianInv(jacobian_inv_cell, jacobian_cell);
427 CellTools<double>::setJacobianDet(jacobian_det_cell, jacobian_cell);
428
429 // compute weighted measure
430 FunctionSpaceTools::computeCellMeasure<double>(weighted_measure_cell, jacobian_det_cell, cub_weights_cell);
431
433 // Computing mass matrices:
434 // tabulate values of basis functions at (reference) cubature points
435 basis->getValues(value_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_VALUE);
436
437 // transform values of basis functions
438 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_cell,
439 value_of_basis_at_cub_points_cell);
440
441 // multiply with weighted measure
442 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_cell,
443 weighted_measure_cell,
444 transformed_value_of_basis_at_cub_points_cell);
445
446 // compute mass matrices
447 FunctionSpaceTools::integrate<double>(fe_matrix,
448 transformed_value_of_basis_at_cub_points_cell,
449 weighted_transformed_value_of_basis_at_cub_points_cell,
450 COMP_BLAS);
452
454 // Computing stiffness matrices:
455 // tabulate gradients of basis functions at (reference) cubature points
456 basis->getValues(grad_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_GRAD);
457
458 // transform gradients of basis functions
459 FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points_cell,
460 jacobian_inv_cell,
461 grad_of_basis_at_cub_points_cell);
462
463 // multiply with weighted measure
464 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points_cell,
465 weighted_measure_cell,
466 transformed_grad_of_basis_at_cub_points_cell);
467
468 // compute stiffness matrices and sum into fe_matrix
469 FunctionSpaceTools::integrate<double>(fe_matrix,
470 transformed_grad_of_basis_at_cub_points_cell,
471 weighted_transformed_grad_of_basis_at_cub_points_cell,
472 COMP_BLAS,
473 true);
475
477 // Computing RHS contributions:
478 // map cell (reference) cubature points to physical space
479 CellTools<double>::mapToPhysicalFrame(cub_points_cell_physical, cub_points_cell, cell_nodes[pcell], cell);
480
481 // evaluate rhs function
482 rhsFunc(rhs_at_cub_points_cell_physical, cub_points_cell_physical, x_order, y_order, z_order);
483
484 // compute rhs
485 FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
486 rhs_at_cub_points_cell_physical,
487 weighted_transformed_value_of_basis_at_cub_points_cell,
488 COMP_BLAS);
489
490 // compute neumann b.c. contributions and adjust rhs
491 sideCub->getCubature(cub_points_side, cub_weights_side);
492 for (unsigned i=0; i<numSides; i++) {
493 // compute geometric cell information
494 CellTools<double>::mapToReferenceSubcell(cub_points_side_refcell, cub_points_side, sideDim, (int)i, cell);
495 CellTools<double>::setJacobian(jacobian_side_refcell, cub_points_side_refcell, cell_nodes[pcell], cell);
496 CellTools<double>::setJacobianDet(jacobian_det_side_refcell, jacobian_side_refcell);
497
498 // compute weighted face measure
499 FunctionSpaceTools::computeFaceMeasure<double>(weighted_measure_side_refcell,
500 jacobian_side_refcell,
501 cub_weights_side,
502 i,
503 cell);
504
505 // tabulate values of basis functions at side cubature points, in the reference parent cell domain
506 basis->getValues(value_of_basis_at_cub_points_side_refcell, cub_points_side_refcell, OPERATOR_VALUE);
507 // transform
508 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_side_refcell,
509 value_of_basis_at_cub_points_side_refcell);
510
511 // multiply with weighted measure
512 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_side_refcell,
513 weighted_measure_side_refcell,
514 transformed_value_of_basis_at_cub_points_side_refcell);
515
516 // compute Neumann data
517 // map side cubature points in reference parent cell domain to physical space
518 CellTools<double>::mapToPhysicalFrame(cub_points_side_physical, cub_points_side_refcell, cell_nodes[pcell], cell);
519 // now compute data
520 neumann(neumann_data_at_cub_points_side_physical, cub_points_side_physical, jacobian_side_refcell,
521 cell, (int)i, x_order, y_order, z_order);
522
523 FunctionSpaceTools::integrate<double>(neumann_fields_per_side,
524 neumann_data_at_cub_points_side_physical,
525 weighted_transformed_value_of_basis_at_cub_points_side_refcell,
526 COMP_BLAS);
527
528 // adjust RHS
529 RealSpaceTools<double>::add(rhs_and_soln_vector, neumann_fields_per_side);;
530 }
532
534 // Solution of linear system:
535 int info = 0;
536 Teuchos::LAPACK<int, double> solver;
537 solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
539
541 // Building interpolant:
542 // evaluate basis at interpolation points
543 basis->getValues(value_of_basis_at_interp_points_ref, interp_points_ref, OPERATOR_VALUE);
544 // transform values of basis functions
545 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points_ref,
546 value_of_basis_at_interp_points_ref);
547 FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points_ref);
549
550 /******************* END COMPUTATION ***********************/
551
552 RealSpaceTools<double>::subtract(interpolant, exact_solution);
553
554 *outStream << "\nRelative norm-2 error between exact solution polynomial of order ("
555 << x_order << ", " << y_order << ", " << z_order
556 << ") and finite element interpolant of order " << basis_order << ": "
557 << RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
558 RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) << "\n";
559
560 if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
561 RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) > zero) {
562 *outStream << "\n\nPatch test failed for solution polynomial order ("
563 << x_order << ", " << y_order << ", " << z_order << ") and basis order " << basis_order << "\n\n";
564 errorFlag++;
565 }
566 } // end for z_order
567 } // end for y_order
568 } // end for x_order
569 } // end for pcell
570
571 }
572 // Catch unexpected errors
573 catch (const std::logic_error & err) {
574 *outStream << err.what() << "\n\n";
575 errorFlag = -1000;
576 };
577
578 if (errorFlag != 0)
579 std::cout << "End Result: TEST FAILED\n";
580 else
581 std::cout << "End Result: TEST PASSED\n";
582
583 // reset format state of std::cout
584 std::cout.copyfmt(oldFormatState);
585
586 return errorFlag;
587}
void rhsFunc(FieldContainer< double > &, const FieldContainer< double > &, int, int, int)
right-hand side function
Definition: test_02.cpp:76
void u_exact(FieldContainer< double > &, const FieldContainer< double > &, int, int, int)
exact solution
Definition: test_02.cpp:181
void neumann(FieldContainer< double > &, const FieldContainer< double > &, const FieldContainer< double > &, const shards::CellTopology &, int, int, int, int)
neumann boundary conditions
Definition: test_02.cpp:125
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_HEX_Cn_FEM class.
Header file for utility class to provide point tools, such as barycentric coordinates,...
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 2 on Hexahedron cell.
A stateless class for operations on cell data. Provides methods for:
static void mapToReferenceSubcell(ArraySubcellPoint &refSubcellPoints, const ArrayParamPoint &paramPoints, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
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.
static void getPhysicalSideNormals(ArraySideNormal &sideNormals, const ArrayJac &worksetJacobians, const int &worksetSideOrd, const shards::CellTopology &parentCell)
Computes non-normalized normal vectors to physical sides in a side workset .
static void setJacobianDet(ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F.
static void setJacobianInv(ArrayJacInv &jacobianInv, const ArrayJac &jacobian)
Computes the inverse of the Jacobian matrix DF of the reference-to-physical frame map F.
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.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
int dimension(const int whichDim) const
Returns the specified dimension.
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...
Implementation of basic linear algebra functionality in Euclidean space.
static void subtract(Scalar *diffArray, const Scalar *inArray1, const Scalar *inArray2, const int size)
Subtracts contiguous data inArray2 from inArray1 of size size: diffArray = inArray1 - inArray2.
static Scalar vectorNorm(const Scalar *inVec, const size_t dim, const ENorm normType)
Computes norm (1, 2, infinity) of the vector inVec of size dim.
static void add(Scalar *sumArray, const Scalar *inArray1, const Scalar *inArray2, const int size)
Adds contiguous data inArray1 and inArray2 of size size: sumArray = inArray1 + inArray2.