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 
60 int 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
119  controlVolSideCub = Teuchos::rcp(new Intrepid::CubatureControlVolumeSide<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo));
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 (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
315  controlVolSideCub = Teuchos::rcp(new Intrepid::CubatureControlVolumeSide<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo));
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 (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
514  controlVolSideCub = Teuchos::rcp(new Intrepid::CubatureControlVolumeSide<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo));
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 (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
723  controlVolSideCub = Teuchos::rcp(new Intrepid::CubatureControlVolumeSide<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo));
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 (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 (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 (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
1016  controlVolCub = Teuchos::rcp(new Intrepid::CubatureControlVolume<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo));
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 (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
1064  controlVolCub = Teuchos::rcp(new Intrepid::CubatureControlVolume<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo));
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 (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
1113  controlVolCub = Teuchos::rcp(new Intrepid::CubatureControlVolume<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo));
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 (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
1217  controlVolCub = Teuchos::rcp(new Intrepid::CubatureControlVolume<double,Intrepid::FieldContainer<double>,Intrepid::FieldContainer<double> >(cellTopo));
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 (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 }
Defines cubature (integration) rules over Neumann boundaries for control volume method.
Intrepid utilities.
Header file for the Intrepid::CubatureControlVolume class.
Header file for the Intrepid::CubatureControlVolumeSide class.
Header file for the Intrepid::CubatureControlVolumeBoundary class.
Defines cubature (integration) rules over control volumes.
Defines cubature (integration) rules over control volumes.