Intrepid
test_25.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
44
54#include "Intrepid_Utils.hpp"
55#include "Teuchos_oblackholestream.hpp"
56#include "Teuchos_RCP.hpp"
57#include "Teuchos_GlobalMPISession.hpp"
58
59
60int main(int argc, char *argv[]) {
61
62 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
63
64 // This little trick lets us print to std::cout only if
65 // a (dummy) command-line argument is provided.
66 int iprint = argc - 1;
67 Teuchos::RCP<std::ostream> outStream;
68 Teuchos::oblackholestream bhs; // outputs nothing
69 if (iprint > 0)
70 outStream = Teuchos::rcp(&std::cout, false);
71 else
72 outStream = Teuchos::rcp(&bhs, false);
73
74 // Save the format state of the original std::cout.
75 Teuchos::oblackholestream oldFormatState;
76 oldFormatState.copyfmt(std::cout);
77
78 *outStream \
79 << "===============================================================================\n" \
80 << "| |\n" \
81 << "| Unit Test (CubatureControlVolume) |\n" \
82 << "| (CubatureControlVolumeSide) |\n" \
83 << "| (CubatureControlVolumeBoundary) |\n" \
84 << "| |\n" \
85 << "| 1) Correctness of cubature points and weights |\n"\
86 << "| 2) Comparison of sub-control volume weights and primary cell volume |\n"\
87 << "| 3) Control volume integration |\n"\
88 << "| |\n" \
89 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
90 << "| Denis Ridzal (dridzal@sandia.gov) or |\n" \
91 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
92 << "| |\n" \
93 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
94 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
95 << "| |\n" \
96 << "===============================================================================\n"\
97 << "| TEST 1: correctness of cubature points and weights |\n"\
98 << "===============================================================================\n";
99
100 int errorFlag = 0;
101
102 Teuchos::RCP<shards::CellTopology> cellTopo;
103 Teuchos::RCP<Intrepid::Cubature<double,Intrepid::FieldContainer<double> > > controlVolCub;
104 Teuchos::RCP<Intrepid::Cubature<double,Intrepid::FieldContainer<double> > > controlVolSideCub;
105 Teuchos::RCP<Intrepid::Cubature<double,Intrepid::FieldContainer<double> > > controlVolBCCub;
106
107
108 // quadrilateral primary cells
109 try {
110
111 // set cell topology
112 cellTopo = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<> >()));
113
114 // define control volume cubature rule
116 int numPoints = controlVolCub->getNumPoints();
117
118 // define control volume side cubature rule
120 int numSidePoints = controlVolSideCub->getNumPoints();
121
122 // define control volume boundary cubature rule
123 int iside = 2; // boundary side of primary cell
124 controlVolBCCub = Teuchos::rcp(new Intrepid::CubatureControlVolumeBoundary<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo,iside));
125 int numBCPoints = controlVolBCCub->getNumPoints();
126
127 // define quad coordinates
128 Intrepid::FieldContainer<double> cellCoords(2,4,2);
129 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0;
130 cellCoords(0,1,0) = 0.5; cellCoords(0,1,1) = 0.0;
131 cellCoords(0,2,0) = 0.5; cellCoords(0,2,1) = 1.0;
132 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 1.0;
133 cellCoords(1,0,0) = 0.5; cellCoords(1,0,1) = 0.0;
134 cellCoords(1,1,0) = 1.0; cellCoords(1,1,1) = 0.0;
135 cellCoords(1,2,0) = 1.0; cellCoords(1,2,1) = 1.0;
136 cellCoords(1,3,0) = 0.5; cellCoords(1,3,1) = 1.0;
137
138 // points
139 double exactCubPoints[] = {
140 0.125, 0.25, 0.375, 0.25,
141 0.375, 0.75, 0.125, 0.75,
142 0.625, 0.25, 0.875, 0.25,
143 0.875, 0.75, 0.625, 0.75
144 };
145
146 // weights
147 double exactCubWeights[] = {
148 0.125, 0.125, 0.125, 0.125,
149 0.125, 0.125, 0.125, 0.125
150 };
151
152 // side points
153 double exactSideCubPoints[] = {
154 0.25, 0.25, 0.375, 0.5,
155 0.25, 0.75, 0.125, 0.5,
156 0.75, 0.25, 0.875, 0.5,
157 0.75, 0.75, 0.625, 0.5
158 };
159
160 // side weights (these are weighted normals!)
161 double exactSideCubWeights[] = {
162 0.5, 0.0, 0.0, 0.25,
163 -0.5, 0.0, 0.0,-0.25,
164 0.5, 0.0, 0.0, 0.25,
165 -0.5, 0.0, 0.0,-0.25
166 };
167
168 // boundary points
169 double exactBCCubPoints[] = {
170 0.375, 1.0, 0.125, 1.0,
171 0.875, 1.0, 0.625, 1.0
172 };
173
174 // boundary weights
175 double exactBCCubWeights[] = {
176 0.25, 0.25, 0.25, 0.25
177 };
178
179 int exactPoints = 4;
180 if (numPoints != exactPoints) {
181 errorFlag++;
182 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
183 *outStream << "Number of cubature points: " << numPoints;
184 *outStream << " does not equal correct number: " << exactPoints << "\n";
185 }
186
187 int exactBCPoints = 2;
188 if (numBCPoints != exactBCPoints) {
189 errorFlag++;
190 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
191 *outStream << "Number of boundary cubature points: " << numBCPoints;
192 *outStream << " does not equal correct number: " << exactBCPoints << "\n";
193 }
194
195 // get cubature points and weights for volume integration over control volume
196 int numCells = 2; int spaceDim = 2;
197 Intrepid::FieldContainer<double> cubPoints(numCells,numPoints,spaceDim);
198 Intrepid::FieldContainer<double> cubWeights(numCells,numPoints);
199 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
200
201 // get cubature points and weights for surface integration over control volume
202 // (Note: the weights here are really weighted normals)
203 Intrepid::FieldContainer<double> sideCubPoints(numCells,numSidePoints,spaceDim);
204 Intrepid::FieldContainer<double> sideCubWeights(numCells,numSidePoints,spaceDim);
205 controlVolSideCub->getCubature(sideCubPoints,sideCubWeights,cellCoords);
206
207 // get cubature points and weights for surface integration over Neumann boundary
208 // (Note: this is a boundary of the primary cell)
209 Intrepid::FieldContainer<double> bcCubPoints(numCells,numBCPoints,spaceDim);
210 Intrepid::FieldContainer<double> bcCubWeights(numCells,numBCPoints);
211 controlVolBCCub->getCubature(bcCubPoints,bcCubWeights,cellCoords);
212
213 int countp = 0;
214 int countw = 0;
215 int countbcp = 0;
216 int countbcw = 0;
217 for (int i=0; i<numCells; i++) {
218 for (int j=0; j<numPoints; j++) {
219 for (int k=0; k<spaceDim; k++) {
220
221 // check values of cubature points
222 if (std::abs(cubPoints(i,j,k) - exactCubPoints[countp]) > Intrepid::INTREPID_TOL) {
223 errorFlag++;
224 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
225 *outStream << " At multi-index { ";
226 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
227 *outStream << "} computed point: " << cubPoints(i,j,k)
228 << " but reference value: " << exactCubPoints[countp] << "\n";
229 }
230
231 // check values of side cubature points
232 if (std::abs(sideCubPoints(i,j,k) - exactSideCubPoints[countp]) > Intrepid::INTREPID_TOL) {
233 errorFlag++;
234 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
235
236 // Output the multi-index of the value where the error is:
237 *outStream << " At multi-index { ";
238 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
239 *outStream << "} computed point: " << sideCubPoints(i,j,k)
240 << " but reference value: " << exactSideCubPoints[countp] << "\n";
241 }
242
243 // check values of side cubature weights
244 if (std::abs(sideCubWeights(i,j,k) - exactSideCubWeights[countp]) > Intrepid::INTREPID_TOL) {
245 errorFlag++;
246 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
247
248 // Output the multi-index of the value where the error is:
249 *outStream << " At multi-index { ";
250 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
251 *outStream << "} computed weight: " << sideCubWeights(i,j,k)
252 << " but reference value: " << exactSideCubWeights[countp] << "\n";
253 }
254
255 countp++;
256
257 }
258
259 // check values of cubature weights
260 if (std::abs(cubWeights(i,j) - exactCubWeights[countw]) > Intrepid::INTREPID_TOL) {
261 errorFlag++;
262 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
263 *outStream << " At multi-index { ";
264 *outStream << i << " ";*outStream << j << " ";
265 *outStream << "} computed weight: " << cubWeights(i,j)
266 << " but reference value: " << exactCubWeights[countw] << "\n";
267 }
268 countw++;
269 }
270
271 for (int j=0; j<numBCPoints; j++) {
272 for (int k=0; k<spaceDim; k++) {
273
274 // check values of boundary cubature points
275 if (std::abs(bcCubPoints(i,j,k) - exactBCCubPoints[countbcp]) > Intrepid::INTREPID_TOL) {
276 errorFlag++;
277 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
278 *outStream << " At multi-index { ";
279 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
280 *outStream << "} computed point: " << bcCubPoints(i,j,k)
281 << " but reference value: " << exactBCCubPoints[countbcp] << "\n";
282 }
283 countbcp++;
284 }
285
286 // check values of cubature weights
287 if (std::abs(bcCubWeights(i,j) - exactBCCubWeights[countbcw]) > Intrepid::INTREPID_TOL) {
288 errorFlag++;
289 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
290 *outStream << " At multi-index { ";
291 *outStream << i << " ";*outStream << j << " ";
292 *outStream << "} computed weight: " << bcCubWeights(i,j)
293 << " but reference value: " << exactBCCubWeights[countbcw] << "\n";
294 }
295 countbcw++;
296 }
297 }
298 }
299 catch (const std::exception & err) {
300 *outStream << err.what() << "\n";
301 errorFlag = -1;
302 }
303
304 // triangle primary cells
305 try {
306
307 // set cell topology
308 cellTopo = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Triangle<> >()));
309
310 // define control volume cubature rule
312 int numPoints = controlVolCub->getNumPoints();
313
314 // define control volume side cubature rule
316 int numSidePoints = controlVolSideCub->getNumPoints();
317
318 // define control volume boundary cubature rule
319 int iside = 1; // boundary side of primary cell
320 controlVolBCCub = Teuchos::rcp(new Intrepid::CubatureControlVolumeBoundary<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo,iside));
321 int numBCPoints = controlVolBCCub->getNumPoints();
322
323 // define triangle coordinates
324 Intrepid::FieldContainer<double> cellCoords(2,3,2);
325 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0;
326 cellCoords(0,1,0) = 0.5; cellCoords(0,1,1) = 0.0;
327 cellCoords(0,2,0) = 0.0; cellCoords(0,2,1) = 0.5;
328 cellCoords(1,0,0) = 0.5; cellCoords(1,0,1) = 0.0;
329 cellCoords(1,1,0) = 0.5; cellCoords(1,1,1) = 0.5;
330 cellCoords(1,2,0) = 0.0; cellCoords(1,2,1) = 0.5;
331
332 // points
333 double exactCubPoints[] = {
334 0.1041666666666667, 0.1041666666666667, 0.2916666666666667, 0.1041666666666667,
335 0.1041666666666667, 0.2916666666666667, 0.3958333333333333, 0.2083333333333333,
336 0.3958333333333333, 0.3958333333333333, 0.2083333333333333, 0.3958333333333333
337 };
338
339 // weights
340 double exactCubWeights[] = {
341 0.0416666666666667, 0.0416666666666667, 0.0416666666666667,
342 0.0416666666666667, 0.0416666666666667, 0.0416666666666667
343 };
344
345 // side points
346 double exactSideCubPoints[] = {
347 0.2083333333333333, 0.0833333333333333, 0.2083333333333333, 0.2083333333333333,
348 0.0833333333333333, 0.2083333333333333, 0.4166666666666667, 0.2916666666666667,
349 0.2916666666666667, 0.4166666666666667, 0.2916666666666667, 0.2916666666666667
350 };
351
352 // side weights (these are weighted normals!)
353 double exactSideCubWeights[] = {
354 0.1666666666666667, 0.0833333333333333,-0.0833333333333333, 0.0833333333333333,
355 -0.0833333333333333,-0.1666666666666667, 0.0833333333333333, 0.1666666666666667,
356 -0.1666666666666667,-0.0833333333333333, 0.0833333333333333,-0.0833333333333333
357 };
358
359 // boundary points
360 double exactBCCubPoints[] = {
361 0.375, 0.125, 0.125, 0.375,
362 0.375, 0.5, 0.125, 0.5
363 };
364
365 // boundary weights
366 double exactBCCubWeights[] = {
367 0.353553390593274, 0.353553390593274, 0.25, 0.25
368 };
369
370 // same number of points for volume and side in 2-D
371 int exactPoints = 3;
372 if (numPoints != exactPoints) {
373 errorFlag++;
374 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
375 *outStream << "Number of cubature points: " << numPoints;
376 *outStream << " does not equal correct number: " << exactPoints << "\n";
377 }
378 if (numSidePoints != exactPoints) {
379 errorFlag++;
380 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
381 *outStream << "Number of side cubature points: " << numSidePoints;
382 *outStream << " does not equal correct number: " << exactPoints << "\n";
383 }
384
385 int exactBCPoints = 2;
386 if (numBCPoints != exactBCPoints) {
387 errorFlag++;
388 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
389 *outStream << "Number of boundary cubature points: " << numBCPoints;
390 *outStream << " does not equal correct number: " << exactBCPoints << "\n";
391 }
392
393 // get cubature points and weights for volume integration over control volume
394 int numCells = 2; int spaceDim = 2;
395 Intrepid::FieldContainer<double> cubPoints(numCells,numPoints,spaceDim);
396 Intrepid::FieldContainer<double> cubWeights(numCells,numPoints);
397 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
398
399 // get cubature points and weights for surface integration over control volume
400 // (Note: the weights here are really weighted normals)
401 Intrepid::FieldContainer<double> sideCubPoints(numCells,numSidePoints,spaceDim);
402 Intrepid::FieldContainer<double> sideCubWeights(numCells,numSidePoints,spaceDim);
403 controlVolSideCub->getCubature(sideCubPoints,sideCubWeights,cellCoords);
404
405 // get cubature points and weights for surface integration over Neumann boundary
406 // (Note: this is a boundary of the primary cell)
407 Intrepid::FieldContainer<double> bcCubPoints(numCells,numBCPoints,spaceDim);
408 Intrepid::FieldContainer<double> bcCubWeights(numCells,numBCPoints);
409 controlVolBCCub->getCubature(bcCubPoints,bcCubWeights,cellCoords);
410
411 int countp = 0;
412 int countw = 0;
413 int countbcp = 0;
414 int countbcw = 0;
415 for (int i=0; i<numCells; i++) {
416 for (int j=0; j<numPoints; j++) {
417 for (int k=0; k<spaceDim; k++) {
418
419 // check values of cubature points
420 if (std::abs(cubPoints(i,j,k) - exactCubPoints[countp]) > Intrepid::INTREPID_TOL) {
421 errorFlag++;
422 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
423 *outStream << " At multi-index { ";
424 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
425 *outStream << "} computed point: " << cubPoints(i,j,k)
426 << " but reference value: " << exactCubPoints[countp] << "\n";
427 }
428
429 // check values of side cubature points
430 if (std::abs(sideCubPoints(i,j,k) - exactSideCubPoints[countp]) > Intrepid::INTREPID_TOL) {
431 errorFlag++;
432 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
433
434 // Output the multi-index of the value where the error is:
435 *outStream << " At multi-index { ";
436 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
437 *outStream << "} computed point: " << sideCubPoints(i,j,k)
438 << " but reference value: " << exactSideCubPoints[countp] << "\n";
439 }
440
441 // check values of side cubature weights
442 if (std::abs(sideCubWeights(i,j,k) - exactSideCubWeights[countp]) > Intrepid::INTREPID_TOL) {
443 errorFlag++;
444 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
445
446 // Output the multi-index of the value where the error is:
447 *outStream << " At multi-index { ";
448 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
449 *outStream << "} computed weight: " << sideCubWeights(i,j,k)
450 << " but reference value: " << exactSideCubWeights[countp] << "\n";
451 }
452
453 countp++;
454
455 }
456
457 // check values of cubature weights
458 if (std::abs(cubWeights(i,j) - exactCubWeights[countw]) > Intrepid::INTREPID_TOL) {
459 errorFlag++;
460 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
461 *outStream << " At multi-index { ";
462 *outStream << i << " ";*outStream << j << " ";
463 *outStream << "} computed weight: " << cubWeights(i,j)
464 << " but reference value: " << exactCubWeights[countw] << "\n";
465 }
466 countw++;
467 }
468
469 for (int j=0; j<numBCPoints; j++) {
470 for (int k=0; k<spaceDim; k++) {
471
472 // check values of boundary cubature points
473 if (std::abs(bcCubPoints(i,j,k) - exactBCCubPoints[countbcp]) > Intrepid::INTREPID_TOL) {
474 errorFlag++;
475 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
476 *outStream << " At multi-index { ";
477 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
478 *outStream << "} computed point: " << bcCubPoints(i,j,k)
479 << " but reference value: " << exactBCCubPoints[countbcp] << "\n";
480 }
481 countbcp++;
482 }
483
484 // check values of cubature weights
485 if (std::abs(bcCubWeights(i,j) - exactBCCubWeights[countbcw]) > Intrepid::INTREPID_TOL) {
486 errorFlag++;
487 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
488 *outStream << " At multi-index { ";
489 *outStream << i << " ";*outStream << j << " ";
490 *outStream << "} computed weight: " << bcCubWeights(i,j)
491 << " but reference value: " << bcCubWeights(i,j) - exactBCCubWeights[countbcw] << "\n";
492 }
493 countbcw++;
494 }
495 }
496
497 }
498 catch (const std::exception & err) {
499 *outStream << err.what() << "\n";
500 errorFlag = -1;
501 }
502
503 // tetrahedral primary cells
504 try {
505
506 // set cell topology
507 cellTopo = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Tetrahedron<> >()));
508
509 // define control volume cubature rule
511 int numPoints = controlVolCub->getNumPoints();
512
513 // define control volume side cubature rule
515 int numSidePoints = controlVolSideCub->getNumPoints();
516
517 // define control volume boundary cubature rule
518 int iside = 0; // boundary side of primary cell
519 controlVolBCCub = Teuchos::rcp(new Intrepid::CubatureControlVolumeBoundary<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo,iside));
520 int numBCPoints = controlVolBCCub->getNumPoints();
521
522 // define tetrahedron coordinates
523 Intrepid::FieldContainer<double> cellCoords(1,4,3);
524 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0; cellCoords(0,0,2) = 0.0;
525 cellCoords(0,1,0) = 1.0; cellCoords(0,1,1) = 0.0; cellCoords(0,1,2) = 0.0;
526 cellCoords(0,2,0) = 0.0; cellCoords(0,2,1) = 1.0; cellCoords(0,2,2) = 0.0;
527 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 0.0; cellCoords(0,3,2) = 1.0;
528
529 // points
530 double exactCubPoints[] = {
531 0.17708333333333333, 0.17708333333333333, 0.17708333333333333,
532 0.46875000000000000, 0.17708333333333333, 0.17708333333333333,
533 0.17708333333333333, 0.46875000000000000, 0.17708333333333333,
534 0.17708333333333333, 0.17708333333333333, 0.46875000000000000
535 };
536
537 // weights
538 double exactCubWeights[] = {
539 0.0416666666666667, 0.0416666666666667,
540 0.0416666666666667, 0.0416666666666667
541 };
542
543 // side points
544 double exactSideCubPoints[] = {
545 0.3541666666666667, 0.1458333333333333, 0.1458333333333333,
546 0.3541666666666667, 0.3541666666666667, 0.1458333333333333,
547 0.1458333333333333, 0.3541666666666667, 0.1458333333333333,
548 0.1458333333333333, 0.1458333333333333, 0.3541666666666667,
549 0.3541666666666667, 0.1458333333333333, 0.3541666666666667,
550 0.1458333333333333, 0.3541666666666667, 0.3541666666666667
551 };
552
553 // side weights (these are weighted normals!)
554 double exactSideCubWeights[] = {
555 0.0833333333333333, 0.0416666666666667, 0.041666666666667,
556 -0.0416666666666667, 0.0416666666666667, 0.000000000000000,
557 -0.0416666666666667,-0.0833333333333333,-0.041666666666667,
558 0.0416666666666667, 0.0416666666666667, 0.083333333333333,
559 -0.0416666666666667, 0.0000000000000000, 0.041666666666667,
560 0.0000000000000000,-0.0416666666666667, 0.041666666666667
561 };
562
563 // boundary points
564 double exactBCCubPoints[] = {
565 0.208333333333333, 0.00, 0.208333333333333,
566 0.583333333333333, 0.00, 0.208333333333333,
567 0.208333333333333, 0.00, 0.583333333333333,
568 };
569
570 // boundary weights
571 double exactBCCubWeights[] = {
572 0.166666666666667, 0.166666666666667, 0.166666666666667
573 };
574
575 // volume points
576 int exactPoints = 4;
577 if (numPoints != exactPoints) {
578 errorFlag++;
579 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
580 *outStream << "Number of cubature points: " << numPoints;
581 *outStream << " does not equal correct number: " << exactPoints << "\n";
582 }
583 // side points
584 exactPoints = 6;
585 if (numSidePoints != exactPoints) {
586 errorFlag++;
587 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
588 *outStream << "Number of side cubature points: " << numSidePoints;
589 *outStream << " does not equal correct number: " << exactPoints << "\n";
590 }
591 // boundary points
592 int exactBCPoints = 3;
593 if (numBCPoints != exactBCPoints) {
594 errorFlag++;
595 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
596 *outStream << "Number of boundary cubature points: " << numBCPoints;
597 *outStream << " does not equal correct number: " << exactBCPoints << "\n";
598 }
599
600 // get cubature points and weights for volume integration over control volume
601 int numCells = 1; int spaceDim = 3;
602 Intrepid::FieldContainer<double> cubPoints(numCells,numPoints,spaceDim);
603 Intrepid::FieldContainer<double> cubWeights(numCells,numPoints);
604 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
605
606 // get cubature points and weights for surface integration over control volume
607 // (Note: the weights here are really weighted normals)
608 Intrepid::FieldContainer<double> sideCubPoints(numCells,numSidePoints,spaceDim);
609 Intrepid::FieldContainer<double> sideCubWeights(numCells,numSidePoints,spaceDim);
610 controlVolSideCub->getCubature(sideCubPoints,sideCubWeights,cellCoords);
611
612 // get cubature points and weights for surface integration over Neumann boundary
613 // (Note: this is a boundary of the primary cell)
614 Intrepid::FieldContainer<double> bcCubPoints(numCells,numBCPoints,spaceDim);
615 Intrepid::FieldContainer<double> bcCubWeights(numCells,numBCPoints);
616 controlVolBCCub->getCubature(bcCubPoints,bcCubWeights,cellCoords);
617
618 int countp = 0;
619 int countw = 0;
620 int countbcp = 0;
621 int countbcw = 0;
622 for (int i=0; i<numCells; i++) {
623 for (int j=0; j<numPoints; j++) {
624 for (int k=0; k<spaceDim; k++) {
625
626 // check values of cubature points
627 if (std::abs(cubPoints(i,j,k) - exactCubPoints[countp]) > Intrepid::INTREPID_TOL) {
628 errorFlag++;
629 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
630 *outStream << " At multi-index { ";
631 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
632 *outStream << "} computed point: " << cubPoints(i,j,k)
633 << " but reference value: " << exactCubPoints[countp] << "\n";
634 }
635
636 // check values of side cubature points
637 if (std::abs(sideCubPoints(i,j,k) - exactSideCubPoints[countp]) > Intrepid::INTREPID_TOL) {
638 errorFlag++;
639 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
640
641 // Output the multi-index of the value where the error is:
642 *outStream << " At multi-index { ";
643 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
644 *outStream << "} computed point: " << sideCubPoints(i,j,k)
645 << " but reference value: " << exactSideCubPoints[countp] << "\n";
646 }
647
648 // check values of side cubature weights
649 if (std::abs(sideCubWeights(i,j,k) - exactSideCubWeights[countp]) > Intrepid::INTREPID_TOL) {
650 errorFlag++;
651 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
652
653 // Output the multi-index of the value where the error is:
654 *outStream << " At multi-index { ";
655 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
656 *outStream << "} computed weight: " << sideCubWeights(i,j,k)
657 << " but reference value: " << exactSideCubWeights[countp] << "\n";
658 }
659
660 countp++;
661
662 }
663
664 // check values of cubature weights
665 if (std::abs(cubWeights(i,j) - exactCubWeights[countw]) > Intrepid::INTREPID_TOL) {
666 errorFlag++;
667 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
668 *outStream << " At multi-index { ";
669 *outStream << i << " ";*outStream << j << " ";
670 *outStream << "} computed weight: " << cubWeights(i,j)
671 << " but reference value: " << exactCubWeights[countw] << "\n";
672 *outStream << cubWeights(i,j)-exactCubWeights[countp] << "\n";
673 }
674 countw++;
675 }
676
677 for (int j=0; j<numBCPoints; j++) {
678 for (int k=0; k<spaceDim; k++) {
679
680 // check values of boundary cubature points
681 if (std::abs(bcCubPoints(i,j,k) - exactBCCubPoints[countbcp]) > Intrepid::INTREPID_TOL) {
682 errorFlag++;
683 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
684 *outStream << " At multi-index { ";
685 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
686 *outStream << "} computed point: " << bcCubPoints(i,j,k)
687 << " but reference value: " << exactBCCubPoints[countbcp] << "\n";
688 }
689 countbcp++;
690 }
691
692 // check values of cubature weights
693 if (std::abs(bcCubWeights(i,j) - exactBCCubWeights[countbcw]) > Intrepid::INTREPID_TOL) {
694 errorFlag++;
695 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
696 *outStream << " At multi-index { ";
697 *outStream << i << " ";*outStream << j << " ";
698 *outStream << "} computed weight: " << bcCubWeights(i,j)
699 << " but reference value: " << exactBCCubWeights[countbcw] << "\n";
700 }
701 countbcw++;
702 }
703 }
704
705
706 }
707 catch (const std::exception & err) {
708 *outStream << err.what() << "\n";
709 errorFlag = -1;
710 }
711
712 // hexahedral primary cells
713 try {
714
715 // set cell topology
716 cellTopo = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Hexahedron<> >()));
717
718 // define control volume cubature rule
720 int numPoints = controlVolCub->getNumPoints();
721
722 // define control volume side cubature rule
724 int numSidePoints = controlVolSideCub->getNumPoints();
725
726 // define control volume boundary cubature rule
727 int iside = 0; // boundary side of primary cell
728 controlVolBCCub = Teuchos::rcp(new Intrepid::CubatureControlVolumeBoundary<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo,iside));
729 int numBCPoints = controlVolBCCub->getNumPoints();
730
731 // define hexahedron coordinates
732 Intrepid::FieldContainer<double> cellCoords(1,8,3);
733 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0; cellCoords(0,0,2) = 0.0;
734 cellCoords(0,1,0) = 2.0; cellCoords(0,1,1) = 0.0; cellCoords(0,1,2) = 0.0;
735 cellCoords(0,2,0) = 2.0; cellCoords(0,2,1) = 1.5; cellCoords(0,2,2) = 0.0;
736 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 1.5; cellCoords(0,3,2) = 0.0;
737 cellCoords(0,4,0) = 0.0; cellCoords(0,4,1) = 0.0; cellCoords(0,4,2) = 1.0;
738 cellCoords(0,5,0) = 2.0; cellCoords(0,5,1) = 0.0; cellCoords(0,5,2) = 1.0;
739 cellCoords(0,6,0) = 2.0; cellCoords(0,6,1) = 1.5; cellCoords(0,6,2) = 1.0;
740 cellCoords(0,7,0) = 0.0; cellCoords(0,7,1) = 1.5; cellCoords(0,7,2) = 1.0;
741
742 // points
743 double exactCubPoints[] = {
744 0.5, 0.375, 0.25, 1.5, 0.375, 0.25,
745 1.5, 1.125, 0.25, 0.5, 1.125, 0.25,
746 0.5, 0.375, 0.75, 1.5, 0.375, 0.75,
747 1.5, 1.125, 0.75, 0.5, 1.125, 0.75
748 };
749
750 // weights
751 double exactCubWeights[] = {
752 0.375, 0.375, 0.375, 0.375,
753 0.375, 0.375, 0.375, 0.375
754 };
755
756 // side points
757 double exactSideCubPoints[] = {
758 1.0, 0.375, 0.25, 1.5, 0.750, 0.25,
759 1.0, 1.125, 0.25, 0.5, 0.750, 0.25,
760 1.0, 0.375, 0.75, 1.5, 0.750, 0.75,
761 1.0, 1.125, 0.75, 0.5, 0.750, 0.75,
762 0.5, 0.375, 0.50, 1.5, 0.375, 0.50,
763 1.5, 1.125, 0.50, 0.5, 1.125, 0.50
764 };
765
766 // side weights (these are weighted normals!)
767 double exactSideCubWeights[] = {
768 0.375, 0.00, 0.00, 0.00, 0.50, 0.00,
769 -0.375, 0.00, 0.00, 0.00,-0.50, 0.00,
770 0.375, 0.00, 0.00, 0.00, 0.50, 0.00,
771 -0.375, 0.00, 0.00, 0.00,-0.50, 0.00,
772 0.000, 0.00, 0.75, 0.00, 0.00, 0.75,
773 0.000, 0.00, 0.75, 0.00, 0.00, 0.75,
774 };
775
776 // volume points
777 int exactPoints = 8;
778 if (numPoints != exactPoints) {
779 errorFlag++;
780 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
781 *outStream << "Number of cubature points: " << numPoints;
782 *outStream << " does not equal correct number: " << exactPoints << "\n";
783 }
784 // side points
785 exactPoints = 12;
786 if (numSidePoints != exactPoints) {
787 errorFlag++;
788 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
789 *outStream << "Number of side cubature points: " << numSidePoints;
790 *outStream << " does not equal correct number: " << exactPoints << "\n";
791 }
792
793 // boundary points
794 double exactBCCubPoints[] = {
795 0.5, 0.00, 0.25,
796 1.5, 0.00, 0.25,
797 1.5, 0.00, 0.75,
798 0.5, 0.00, 0.75,
799 };
800
801 // boundary weights
802 double exactBCCubWeights[] = {
803 0.5, 0.5, 0.5, 0.5
804 };
805
806 // get cubature points and weights for volume integration over control volume
807 int numCells = 1; int spaceDim = 3;
808 Intrepid::FieldContainer<double> cubPoints(numCells,numPoints,spaceDim);
809 Intrepid::FieldContainer<double> cubWeights(numCells,numPoints);
810 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
811
812 // get cubature points and weights for surface integration over control volume
813 // (Note: the weights here are really weighted normals)
814 Intrepid::FieldContainer<double> sideCubPoints(numCells,numSidePoints,spaceDim);
815 Intrepid::FieldContainer<double> sideCubWeights(numCells,numSidePoints,spaceDim);
816 controlVolSideCub->getCubature(sideCubPoints,sideCubWeights,cellCoords);
817
818 // get cubature points and weights for surface integration over Neumann boundary
819 // (Note: this is a boundary of the primary cell)
820 Intrepid::FieldContainer<double> bcCubPoints(numCells,numBCPoints,spaceDim);
821 Intrepid::FieldContainer<double> bcCubWeights(numCells,numBCPoints);
822 controlVolBCCub->getCubature(bcCubPoints,bcCubWeights,cellCoords);
823
824 int countp = 0;
825 int countw = 0;
826 int countbcp = 0;
827 int countbcw = 0;
828 for (int i=0; i<numCells; i++) {
829 for (int j=0; j<numPoints; j++) {
830 for (int k=0; k<spaceDim; k++) {
831
832 // check values of cubature points
833 if (std::abs(cubPoints(i,j,k) - exactCubPoints[countp]) > Intrepid::INTREPID_TOL) {
834 errorFlag++;
835 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
836 *outStream << " At multi-index { ";
837 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
838 *outStream << "} computed point: " << cubPoints(i,j,k)
839 << " but reference value: " << exactCubPoints[countp] << "\n";
840 }
841
842 // check values of side cubature points
843 if (std::abs(sideCubPoints(i,j,k) - exactSideCubPoints[countp]) > Intrepid::INTREPID_TOL) {
844 errorFlag++;
845 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
846
847 // Output the multi-index of the value where the error is:
848 *outStream << " At multi-index { ";
849 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
850 *outStream << "} computed point: " << sideCubPoints(i,j,k)
851 << " but reference value: " << exactSideCubPoints[countp] << "\n";
852 }
853
854 // check values of side cubature weights
855 if (std::abs(sideCubWeights(i,j,k) - exactSideCubWeights[countp]) > Intrepid::INTREPID_TOL) {
856 errorFlag++;
857 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
858
859 // Output the multi-index of the value where the error is:
860 *outStream << " At multi-index { ";
861 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
862 *outStream << "} computed weight: " << sideCubWeights(i,j,k)
863 << " but reference value: " << exactSideCubWeights[countp] << "\n";
864 }
865
866 countp++;
867
868 }
869
870 // check values of cubature weights
871 if (std::abs(cubWeights(i,j) - exactCubWeights[countw]) > Intrepid::INTREPID_TOL) {
872 errorFlag++;
873 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
874 *outStream << " At multi-index { ";
875 *outStream << i << " ";*outStream << j << " ";
876 *outStream << "} computed weight: " << cubWeights(i,j)
877 << " but reference value: " << exactCubWeights[countw] << "\n";
878 *outStream << cubWeights(i,j)-exactCubWeights[countp] << "\n";
879 }
880 countw++;
881 }
882
883 for (int j=0; j<numBCPoints; j++) {
884 for (int k=0; k<spaceDim; k++) {
885
886 // check values of boundary cubature points
887 if (std::abs(bcCubPoints(i,j,k) - exactBCCubPoints[countbcp]) > Intrepid::INTREPID_TOL) {
888 errorFlag++;
889 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
890 *outStream << " At multi-index { ";
891 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
892 *outStream << "} computed point: " << bcCubPoints(i,j,k)
893 << " but reference value: " << exactBCCubPoints[countbcp] << "\n";
894 }
895 countbcp++;
896 }
897
898 // check values of cubature weights
899 if (std::abs(bcCubWeights(i,j) - exactBCCubWeights[countbcw]) > Intrepid::INTREPID_TOL) {
900 errorFlag++;
901 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
902 *outStream << " At multi-index { ";
903 *outStream << i << " ";*outStream << j << " ";
904 *outStream << "} computed weight: " << bcCubWeights(i,j)
905 << " but reference value: " << exactBCCubWeights[countbcw] << "\n";
906 }
907 countbcw++;
908 }
909 }
910
911 }
912 catch (const std::exception & err) {
913 *outStream << err.what() << "\n";
914 errorFlag = -1;
915 }
916
917 *outStream \
918 << "===============================================================================\n"\
919 << "| TEST 2: comparison of sub-control volume weights and primary cell volume |\n"\
920 << "===============================================================================\n";
921
922 // quadrilateral primary cells
923 try {
924
925 // set cell topology
926 cellTopo = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<> >()));
927
928 // define control volume cubature rule
930 int numPoints = controlVolCub->getNumPoints();
931
932 // define quad coordinates
933 Intrepid::FieldContainer<double> cellCoords(1,4,2);
934 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0;
935 cellCoords(0,1,0) = 2.4; cellCoords(0,1,1) = 0.0;
936 cellCoords(0,2,0) = 2.4; cellCoords(0,2,1) = 3.1;
937 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 3.1;
938
939 double exact_area = 2.4*3.1;
940
941 // get cubature points and weights
942 int numCells = 1; int spaceDim = 2;
943 Intrepid::FieldContainer<double> cubPoints(numCells,numPoints,spaceDim);
944 Intrepid::FieldContainer<double> cubWeights(numCells,numPoints);
945 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
946
947 // loop over number of points (equals number of sub-control volumes) and check total volume
948 double total_area = 0.0;
949 for (int i=0; i<numPoints; i++) {
950 total_area += cubWeights(0,i);
951 }
952 if (std::abs(total_area - exact_area) > Intrepid::INTREPID_TOL) {
953 errorFlag++;
954 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
955 *outStream << " Sum of sub-control volume areas: ";
956 *outStream << total_area;
957 *outStream << " does not equal quad primary cell area: " << exact_area << "\n";
958 }
959
960 }
961 catch (const std::exception & err) {
962 *outStream << err.what() << "\n";
963 errorFlag = -1;
964 }
965
966 // triangle primary cells
967 try {
968
969 // set cell topology
970 cellTopo = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Triangle<> >()));
971
972 // define control volume cubature rule
974 int numPoints = controlVolCub->getNumPoints();
975
976 // define quad coordinates
977 Intrepid::FieldContainer<double> cellCoords(1,3,2);
978 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0;
979 cellCoords(0,1,0) = 3.6; cellCoords(0,1,1) = 0.0;
980 cellCoords(0,2,0) = 0.0; cellCoords(0,2,1) = 2.8;
981
982 double exact_area = 0.5*3.6*2.8;
983
984 // get cubature points and weights
985 int numCells = 1; int spaceDim = 2;
986 Intrepid::FieldContainer<double> cubPoints(numCells,numPoints,spaceDim);
987 Intrepid::FieldContainer<double> cubWeights(numCells,numPoints);
988 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
989
990 // loop over number of points (equals number of sub-control volumes) and check total volume
991 double total_area = 0.0;
992 for (int i=0; i<numPoints; i++) {
993 total_area += cubWeights(0,i);
994 }
995 if (std::abs(total_area - exact_area) > Intrepid::INTREPID_TOL) {
996 errorFlag++;
997 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
998 *outStream << " Sum of sub-control volume areas: ";
999 *outStream << total_area;
1000 *outStream << " does not equal triangle primary cell area: " << exact_area << "\n";
1001 }
1002
1003 }
1004 catch (const std::exception & err) {
1005 *outStream << err.what() << "\n";
1006 errorFlag = -1;
1007 }
1008
1009 // hexahedral primary cells
1010 try {
1011
1012 // set cell topology
1013 cellTopo = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Hexahedron<> >()));
1014
1015 // define control volume cubature rule
1017 int numPoints = controlVolCub->getNumPoints();
1018
1019 // define hexahedron coordinates
1020 Intrepid::FieldContainer<double> cellCoords(1,8,3);
1021 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0; cellCoords(0,0,2) = 0.0;
1022 cellCoords(0,1,0) = 2.4; cellCoords(0,1,1) = 0.0; cellCoords(0,1,2) = 0.0;
1023 cellCoords(0,2,0) = 2.4; cellCoords(0,2,1) = 3.1; cellCoords(0,2,2) = 0.0;
1024 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 3.1; cellCoords(0,3,2) = 0.0;
1025 cellCoords(0,4,0) = 0.0; cellCoords(0,4,1) = 0.0; cellCoords(0,4,2) = 1.7;
1026 cellCoords(0,5,0) = 2.4; cellCoords(0,5,1) = 0.0; cellCoords(0,5,2) = 1.7;
1027 cellCoords(0,6,0) = 2.4; cellCoords(0,6,1) = 3.1; cellCoords(0,6,2) = 1.7;
1028 cellCoords(0,7,0) = 0.0; cellCoords(0,7,1) = 3.1; cellCoords(0,7,2) = 1.7;
1029
1030 double exact_vol = 2.4*3.1*1.7;
1031
1032 // get cubature points and weights
1033 int numCells = 1; int spaceDim = 3;
1034 Intrepid::FieldContainer<double> cubPoints(numCells,numPoints,spaceDim);
1035 Intrepid::FieldContainer<double> cubWeights(numCells,numPoints);
1036 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
1037
1038 // loop over number of points (equals number of sub-control volumes) and check total volume
1039 double total_vol = 0.0;
1040 for (int i=0; i<numPoints; i++) {
1041 total_vol += cubWeights(0,i);
1042 }
1043 if (std::abs(total_vol - exact_vol) > Intrepid::INTREPID_TOL) {
1044 errorFlag++;
1045 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
1046 *outStream << " Sum of sub-control volumes: ";
1047 *outStream << total_vol;
1048 *outStream << " does not equal hexahedron primary cell volume: " << exact_vol << "\n";
1049 }
1050
1051 }
1052 catch (const std::exception & err) {
1053 *outStream << err.what() << "\n";
1054 errorFlag = -1;
1055 }
1056
1057 // tetrahedral primary cells
1058 try {
1059
1060 // set cell topology
1061 cellTopo = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Tetrahedron<> >()));
1062
1063 // define control volume cubature rule
1065 int numPoints = controlVolCub->getNumPoints();
1066
1067 // define tetrahedron coordinates
1068 Intrepid::FieldContainer<double> cellCoords(1,4,3);
1069 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0; cellCoords(0,0,2) = 0.0;
1070 cellCoords(0,1,0) = 3.6; cellCoords(0,1,1) = 0.0; cellCoords(0,1,2) = 0.0;
1071 cellCoords(0,2,0) = 0.0; cellCoords(0,2,1) = 2.8; cellCoords(0,2,2) = 0.0;
1072 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 2.8; cellCoords(0,3,2) = 1.7;
1073
1074 double exact_vol = (0.5*3.6*2.8)*1.7/3.0;
1075
1076 // get cubature points and weights
1077 int numCells = 1; int spaceDim = 3;
1078 Intrepid::FieldContainer<double> cubPoints(numCells,numPoints,spaceDim);
1079 Intrepid::FieldContainer<double> cubWeights(numCells,numPoints);
1080 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
1081
1082 // loop over number of points (equals number of sub-control volumes) and check total volume
1083 double total_vol = 0.0;
1084 for (int i=0; i<numPoints; i++) {
1085 total_vol += cubWeights(0,i);
1086 }
1087 if (std::abs(total_vol - exact_vol) > Intrepid::INTREPID_TOL) {
1088 errorFlag++;
1089 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
1090 *outStream << " Sum of sub-control volumes: ";
1091 *outStream << total_vol;
1092 *outStream << " does not equal tetrahedron primary cell volume: " << exact_vol << "\n";
1093 }
1094
1095 }
1096 catch (const std::exception & err) {
1097 *outStream << err.what() << "\n";
1098 errorFlag = -1;
1099 }
1100
1101 *outStream \
1102 << "===============================================================================\n"\
1103 << "| TEST 3: control volume integration |\n"\
1104 << "===============================================================================\n";
1105
1106 // quadrilateral primary cells
1107 try {
1108
1109 // set cell topology
1110 cellTopo = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<> >()));
1111
1112 // define control volume cubature rule
1114 int numPoints = controlVolCub->getNumPoints();
1115
1116 // define control volume cubature rule
1117 controlVolSideCub = Teuchos::rcp(new Intrepid::CubatureControlVolumeSide<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo));
1118 int numSidePoints = controlVolSideCub->getNumPoints();
1119
1120 // define quad coordinates - four cells that define a complete control volume around center node
1121 Intrepid::FieldContainer<double> cellCoords(4,4,2);
1122 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0;
1123 cellCoords(0,1,0) = 0.5; cellCoords(0,1,1) = 0.0;
1124 cellCoords(0,2,0) = 0.41; cellCoords(0,2,1) = 0.58;
1125 cellCoords(0,3,0) = 0.0; cellCoords(0,3,1) = 0.5;
1126 cellCoords(1,0,0) = 0.5; cellCoords(1,0,1) = 0.0;
1127 cellCoords(1,1,0) = 1.0; cellCoords(1,1,1) = 0.0;
1128 cellCoords(1,2,0) = 1.0; cellCoords(1,2,1) = 0.5;
1129 cellCoords(1,3,0) = 0.41; cellCoords(1,3,1) = 0.58;
1130 cellCoords(2,0,0) = 0.0; cellCoords(2,0,1) = 0.5;
1131 cellCoords(2,1,0) = 0.41; cellCoords(2,1,1) = 0.58;
1132 cellCoords(2,2,0) = 0.5; cellCoords(2,2,1) = 1.0;
1133 cellCoords(2,3,0) = 0.0; cellCoords(2,3,1) = 1.0;
1134 cellCoords(3,0,0) = 0.41; cellCoords(3,0,1) = 0.58;
1135 cellCoords(3,1,0) = 1.0; cellCoords(3,1,1) = 0.5;
1136 cellCoords(3,2,0) = 1.0; cellCoords(3,2,1) = 1.0;
1137 cellCoords(3,3,0) = 0.5; cellCoords(3,3,1) = 1.0;
1138
1139 // get cubature points and weights
1140 int numCells = 4; int spaceDim = 2;
1141 Intrepid::FieldContainer<double> cubPoints(numCells,numPoints,spaceDim);
1142 Intrepid::FieldContainer<double> cubWeights(numCells,numPoints);
1143 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
1144
1145 Intrepid::FieldContainer<double> sideCubPoints(numCells,numSidePoints,spaceDim);
1146 Intrepid::FieldContainer<double> sideCubWeights(numCells,numSidePoints,spaceDim);
1147 controlVolSideCub->getCubature(sideCubPoints,sideCubWeights,cellCoords);
1148
1149 // test cubature rule by checking the equality of \int_C \nabla \cdot F dV = \int_dC F \cdot n dS
1150 // using F = (a x,b y), (\nabla \cdot F) = a + b.
1151
1152 // first evaluate F at all control volume side points
1153 double a = 2.1; double b = 1.4;
1154 Intrepid::FieldContainer<double> F(numCells,numSidePoints,spaceDim);
1155 for (int i = 0; i < numCells; i++) {
1156 for (int j = 0; j < numSidePoints; j++) {
1157 F(i,j,0) = a*sideCubPoints(i,j,0);
1158 F(i,j,1) = b*sideCubPoints(i,j,1);
1159 }
1160 }
1161
1162 // gather the correct contributions to the surface integral
1163 double surfaceInt = 0.0;
1164
1165 // contributions from first cell
1166 surfaceInt += - F(0,1,0)*sideCubWeights(0,1,0) - F(0,1,1)*sideCubWeights(0,1,1)
1167 + F(0,2,0)*sideCubWeights(0,2,0) + F(0,2,1)*sideCubWeights(0,2,1);
1168
1169 // contributions from second cell
1170 surfaceInt += - F(1,2,0)*sideCubWeights(1,2,0) - F(1,2,1)*sideCubWeights(1,2,1)
1171 + F(1,3,0)*sideCubWeights(1,3,0) + F(1,3,1)*sideCubWeights(1,3,1);
1172
1173 // contributions from third cell
1174 surfaceInt += - F(2,0,0)*sideCubWeights(2,0,0) - F(2,0,1)*sideCubWeights(2,0,1)
1175 + F(2,1,0)*sideCubWeights(2,1,0) + F(2,1,1)*sideCubWeights(2,1,1);
1176
1177 // contributions from fourth cell
1178 surfaceInt += - F(3,3,0)*sideCubWeights(3,3,0) - F(3,3,1)*sideCubWeights(3,3,1)
1179 + F(3,0,0)*sideCubWeights(3,0,0) + F(3,0,1)*sideCubWeights(3,0,1);
1180
1181 // gather the correct contributions to the volume integral
1182 double volumeInt = 0.0;
1183
1184 // contributions from first cell
1185 volumeInt += (a+b)*cubWeights(0,2);
1186
1187 // contributions from second cell
1188 volumeInt += (a+b)*cubWeights(1,3);
1189
1190 // contributions from third cell
1191 volumeInt += (a+b)*cubWeights(2,1);
1192
1193 // contributions from fourth cell
1194 volumeInt += (a+b)*cubWeights(3,0);
1195
1196 if (std::abs(surfaceInt - volumeInt) > Intrepid::INTREPID_TOL) {
1197 errorFlag++;
1198 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
1199 *outStream << " Integral of (F cdot n) over surface : ";
1200 *outStream << surfaceInt;
1201 *outStream << " does not equal integral of div F over volume: " << volumeInt << "\n";
1202 }
1203
1204 }
1205 catch (const std::exception & err) {
1206 *outStream << err.what() << "\n";
1207 errorFlag = -1;
1208 }
1209
1210 // triangle primary cells
1211 try {
1212
1213 // set cell topology
1214 cellTopo = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Triangle<> >()));
1215
1216 // define control volume cubature rule
1218 int numPoints = controlVolCub->getNumPoints();
1219
1220 // define control volume cubature rule
1221 controlVolSideCub = Teuchos::rcp(new Intrepid::CubatureControlVolumeSide<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo));
1222 int numSidePoints = controlVolSideCub->getNumPoints();
1223
1224 // define triangle coordinates - four cells that define a complete control volume around center node
1225 Intrepid::FieldContainer<double> cellCoords(4,3,2);
1226 cellCoords(0,0,0) = 0.0; cellCoords(0,0,1) = 0.0;
1227 cellCoords(0,1,0) = 1.0; cellCoords(0,1,1) = 0.0;
1228 cellCoords(0,2,0) = 0.41; cellCoords(0,2,1) = 0.58;
1229 cellCoords(1,0,0) = 1.0; cellCoords(1,0,1) = 0.0;
1230 cellCoords(1,1,0) = 1.0; cellCoords(1,1,1) = 1.0;
1231 cellCoords(1,2,0) = 0.41; cellCoords(1,2,1) = 0.58;
1232 cellCoords(2,0,0) = 1.0; cellCoords(2,0,1) = 1.0;
1233 cellCoords(2,1,0) = 0.0; cellCoords(2,1,1) = 1.0;
1234 cellCoords(2,2,0) = 0.41; cellCoords(2,2,1) = 0.58;
1235 cellCoords(3,0,0) = 0.0; cellCoords(3,0,1) = 1.0;
1236 cellCoords(3,1,0) = 0.0; cellCoords(3,1,1) = 0.0;
1237 cellCoords(3,2,0) = 0.41; cellCoords(3,2,1) = 0.58;
1238
1239 // get cubature points and weights
1240 int numCells = 4; int spaceDim = 2;
1241 Intrepid::FieldContainer<double> cubPoints(numCells,numPoints,spaceDim);
1242 Intrepid::FieldContainer<double> cubWeights(numCells,numPoints);
1243 controlVolCub->getCubature(cubPoints,cubWeights,cellCoords);
1244
1245 Intrepid::FieldContainer<double> sideCubPoints(numCells,numSidePoints,spaceDim);
1246 Intrepid::FieldContainer<double> sideCubWeights(numCells,numSidePoints,spaceDim);
1247 controlVolSideCub->getCubature(sideCubPoints,sideCubWeights,cellCoords);
1248
1249 // test cubature rule by checking the equality of \int_C \nabla \cdot F dV = \int_dC F \cdot n dS
1250 // using F = (a x,b y), (\nabla \cdot F) = a + b.
1251
1252 // first evaluate F at all control volume side points
1253 double a = 2.1; double b = 1.4;
1254 Intrepid::FieldContainer<double> F(numCells,numSidePoints,spaceDim);
1255 for (int i = 0; i < numCells; i++) {
1256 for (int j = 0; j < numSidePoints; j++) {
1257 F(i,j,0) = a*sideCubPoints(i,j,0);
1258 F(i,j,1) = b*sideCubPoints(i,j,1);
1259 }
1260 }
1261
1262 // gather the correct contributions to the surface integral
1263 double surfaceInt = 0.0;
1264
1265 // contributions from first cell
1266 surfaceInt += - F(0,1,0)*sideCubWeights(0,1,0) - F(0,1,1)*sideCubWeights(0,1,1)
1267 + F(0,2,0)*sideCubWeights(0,2,0) + F(0,2,1)*sideCubWeights(0,2,1);
1268
1269 // contributions from second cell
1270 surfaceInt += - F(1,1,0)*sideCubWeights(1,1,0) - F(1,1,1)*sideCubWeights(1,1,1)
1271 + F(1,2,0)*sideCubWeights(1,2,0) + F(1,2,1)*sideCubWeights(1,2,1);
1272
1273 // contributions from third cell
1274 surfaceInt += - F(2,1,0)*sideCubWeights(2,1,0) - F(2,1,1)*sideCubWeights(2,1,1)
1275 + F(2,2,0)*sideCubWeights(2,2,0) + F(2,2,1)*sideCubWeights(2,2,1);
1276
1277 // contributions from fourth cell
1278 surfaceInt += - F(3,1,0)*sideCubWeights(3,1,0) - F(3,1,1)*sideCubWeights(3,1,1)
1279 + F(3,2,0)*sideCubWeights(3,2,0) + F(3,2,1)*sideCubWeights(3,2,1);
1280
1281 // gather the correct contributions to the volume integral
1282 double volumeInt = 0.0;
1283
1284 // contributions from first cell
1285 volumeInt += (a+b)*cubWeights(0,2);
1286
1287 // contributions from second cell
1288 volumeInt += (a+b)*cubWeights(1,2);
1289
1290 // contributions from third cell
1291 volumeInt += (a+b)*cubWeights(2,2);
1292
1293 // contributions from fourth cell
1294 volumeInt += (a+b)*cubWeights(3,2);
1295
1296 if (std::abs(surfaceInt - volumeInt) > Intrepid::INTREPID_TOL) {
1297 errorFlag++;
1298 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
1299 *outStream << " Integral of (F cdot n) over surface : ";
1300 *outStream << surfaceInt;
1301 *outStream << " does not equal integral of div F over volume: " << volumeInt << "\n";
1302 }
1303
1304 }
1305 catch (const std::exception & err) {
1306 *outStream << err.what() << "\n";
1307 errorFlag = -1;
1308 }
1309
1310
1311 if (errorFlag != 0)
1312 std::cout << "End Result: TEST FAILED\n";
1313 else
1314 std::cout << "End Result: TEST PASSED\n";
1315
1316 // reset format state of std::cout
1317 std::cout.copyfmt(oldFormatState);
1318
1319 return errorFlag;
1320}
Header file for the Intrepid::CubatureControlVolumeBoundary class.
Header file for the Intrepid::CubatureControlVolumeSide class.
Header file for the Intrepid::CubatureControlVolume class.
static const double INTREPID_TOL
General purpose tolerance in, e.g., internal Newton's method to invert ref to phys maps.
Intrepid utilities.
Defines cubature (integration) rules over Neumann boundaries for control volume method.
Defines cubature (integration) rules over control volumes.
Defines cubature (integration) rules over control volumes.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....