Intrepid
example_01.cpp
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 
92 // Intrepid includes
93 #include "Intrepid_FunctionSpaceTools.hpp"
94 #include "Intrepid_FieldContainer.hpp"
95 #include "Intrepid_CellTools.hpp"
96 #include "Intrepid_ArrayTools.hpp"
97 #include "Intrepid_HCURL_HEX_I1_FEM.hpp"
98 #include "Intrepid_HGRAD_HEX_C1_FEM.hpp"
99 #include "Intrepid_RealSpaceTools.hpp"
100 #include "Intrepid_DefaultCubatureFactory.hpp"
101 #include "Intrepid_Utils.hpp"
102 
103 // Epetra includes
104 #include "Epetra_Time.h"
105 #include "Epetra_Map.h"
106 #include "Epetra_SerialComm.h"
107 #include "Epetra_FECrsMatrix.h"
108 #include "Epetra_FEVector.h"
109 #include "Epetra_Vector.h"
110 
111 // Teuchos includes
112 #include "Teuchos_oblackholestream.hpp"
113 #include "Teuchos_RCP.hpp"
114 #include "Teuchos_BLAS.hpp"
115 
116 // Shards includes
117 #include "Shards_CellTopology.hpp"
118 
119 // EpetraExt includes
120 #include "EpetraExt_RowMatrixOut.h"
121 #include "EpetraExt_MultiVectorOut.h"
122 
123 using namespace std;
124 using namespace Intrepid;
125 
126 // Functions to evaluate exact solution and derivatives
127 int evalu(double & uExact0, double & uExact1, double & uExact2, double & x, double & y, double & z);
128 double evalDivu(double & x, double & y, double & z);
129 int evalCurlu(double & curlu0, double & curlu1, double & curlu2, double & x, double & y, double & z);
130 int evalGradDivu(double & gradDivu0, double & gradDivu1, double & gradDivu2, double & x, double & y, double & z);
131 
132 int main(int argc, char *argv[]) {
133 
134  //Check number of arguments
135  if (argc < 13) {
136  std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
137  std::cout <<"Usage:\n\n";
138  std::cout <<" ./Intrepid_example_Drivers_Example_01.exe NX NY NZ randomMesh mu1 mu2 mu1LX mu1RX mu1LY mu1RY mu1LZ mu1RZ verbose\n\n";
139  std::cout <<" where \n";
140  std::cout <<" int NX - num intervals in x direction (assumed box domain, -1,1) \n";
141  std::cout <<" int NY - num intervals in y direction (assumed box domain, -1,1) \n";
142  std::cout <<" int NZ - num intervals in z direction (assumed box domain, -1,1) \n";
143  std::cout <<" int randomMesh - 1 if mesh randomizer is to be used 0 if not \n";
144  std::cout <<" double mu1 - material property value for region 1 \n";
145  std::cout <<" double mu2 - material property value for region 2 \n";
146  std::cout <<" double mu1LX - left X boundary for region 1 \n";
147  std::cout <<" double mu1RX - right X boundary for region 1 \n";
148  std::cout <<" double mu1LY - left Y boundary for region 1 \n";
149  std::cout <<" double mu1RY - right Y boundary for region 1 \n";
150  std::cout <<" double mu1LZ - bottom Z boundary for region 1 \n";
151  std::cout <<" double mu1RZ - top Z boundary for region 1 \n";
152  std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
153  exit(1);
154  }
155 
156  // This little trick lets us print to std::cout only if
157  // a (dummy) command-line argument is provided.
158  int iprint = argc - 1;
159  Teuchos::RCP<std::ostream> outStream;
160  Teuchos::oblackholestream bhs; // outputs nothing
161  if (iprint > 12)
162  outStream = Teuchos::rcp(&std::cout, false);
163  else
164  outStream = Teuchos::rcp(&bhs, false);
165 
166  // Save the format state of the original std::cout.
167  Teuchos::oblackholestream oldFormatState;
168  oldFormatState.copyfmt(std::cout);
169 
170  *outStream \
171  << "===============================================================================\n" \
172  << "| |\n" \
173  << "| Example: Generate Mass and Stiffness Matrices and Right-Hand Side Vector |\n"
174  << "| for Div-Curl System on Hexahedral Mesh with Curl-Conforming Elements |\n" \
175  << "| |\n" \
176  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
177  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
178  << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
179  << "| |\n" \
180  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
181  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
182  << "| |\n" \
183  << "===============================================================================\n";
184 
185 
186 // ************************************ GET INPUTS **************************************
187 
188  /* In the implementation for discontinuous material properties only the boundaries for
189  region 1, associated with mu1, are input. The remainder of the grid is assumed to use mu2.
190  Note that the material properties are assigned using the undeformed grid. */
191 
192  int NX = atoi(argv[1]); // num intervals in x direction (assumed box domain, -1,1)
193  int NY = atoi(argv[2]); // num intervals in y direction (assumed box domain, -1,1)
194  int NZ = atoi(argv[3]); // num intervals in z direction (assumed box domain, -1,1)
195  int randomMesh = atoi(argv[4]); // 1 if mesh randomizer is to be used 0 if not
196  double mu1 = atof(argv[5]); // material property value for region 1
197  double mu2 = atof(argv[6]); // material property value for region 2
198  double mu1LeftX = atof(argv[7]); // left X boundary for region 1
199  double mu1RightX = atof(argv[8]); // right X boundary for region 1
200  double mu1LeftY = atof(argv[9]); // left Y boundary for region 1
201  double mu1RightY = atof(argv[10]); // right Y boundary for region 1
202  double mu1LeftZ = atof(argv[11]); // left Z boundary for region 1
203  double mu1RightZ = atof(argv[12]); // right Z boundary for region 1
204 
205 // *********************************** CELL TOPOLOGY **********************************
206 
207  // Get cell topology for base hexahedron
208  typedef shards::CellTopology CellTopology;
209  CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() );
210 
211  // Get dimensions
212  int numNodesPerElem = hex_8.getNodeCount();
213  int numEdgesPerElem = hex_8.getEdgeCount();
214  int numFacesPerElem = hex_8.getSideCount();
215  int numNodesPerEdge = 2;
216  int numNodesPerFace = 4;
217  int numEdgesPerFace = 4;
218  int spaceDim = hex_8.getDimension();
219 
220  // Build reference element edge to node map
221  FieldContainer<int> refEdgeToNode(numEdgesPerElem,numNodesPerEdge);
222  for (int i=0; i<numEdgesPerElem; i++){
223  refEdgeToNode(i,0)=hex_8.getNodeMap(1, i, 0);
224  refEdgeToNode(i,1)=hex_8.getNodeMap(1, i, 1);
225  }
226 
227 // *********************************** GENERATE MESH ************************************
228 
229  *outStream << "Generating mesh ... \n\n";
230 
231  *outStream << " NX" << " NY" << " NZ\n";
232  *outStream << std::setw(5) << NX <<
233  std::setw(5) << NY <<
234  std::setw(5) << NZ << "\n\n";
235 
236  // Print mesh information
237  int numElems = NX*NY*NZ;
238  int numNodes = (NX+1)*(NY+1)*(NZ+1);
239  int numEdges = (NX)*(NY + 1)*(NZ + 1) + (NX + 1)*(NY)*(NZ + 1) + (NX + 1)*(NY + 1)*(NZ);
240  int numFaces = (NX)*(NY)*(NZ + 1) + (NX)*(NY + 1)*(NZ) + (NX + 1)*(NY)*(NZ);
241  *outStream << " Number of Elements: " << numElems << " \n";
242  *outStream << " Number of Nodes: " << numNodes << " \n";
243  *outStream << " Number of Edges: " << numEdges << " \n";
244  *outStream << " Number of Faces: " << numFaces << " \n\n";
245 
246  // Cube
247  double leftX = -1.0, rightX = 1.0;
248  double leftY = -1.0, rightY = 1.0;
249  double leftZ = -1.0, rightZ = 1.0;
250 
251  // Mesh spacing
252  double hx = (rightX-leftX)/((double)NX);
253  double hy = (rightY-leftY)/((double)NY);
254  double hz = (rightZ-leftZ)/((double)NZ);
255 
256  // Get nodal coordinates
257  FieldContainer<double> nodeCoord(numNodes, spaceDim);
258  FieldContainer<int> nodeOnBoundary(numNodes);
259  int inode = 0;
260  for (int k=0; k<NZ+1; k++) {
261  for (int j=0; j<NY+1; j++) {
262  for (int i=0; i<NX+1; i++) {
263  nodeCoord(inode,0) = leftX + (double)i*hx;
264  nodeCoord(inode,1) = leftY + (double)j*hy;
265  nodeCoord(inode,2) = leftZ + (double)k*hz;
266  if (k==0 || j==0 || i==0 || k==NZ || j==NY || i==NX){
267  nodeOnBoundary(inode)=1;
268  }
269  inode++;
270  }
271  }
272  }
273 
274  // Element to Node map
275  FieldContainer<int> elemToNode(numElems, numNodesPerElem);
276  int ielem = 0;
277  for (int k=0; k<NZ; k++) {
278  for (int j=0; j<NY; j++) {
279  for (int i=0; i<NX; i++) {
280  elemToNode(ielem,0) = (NY + 1)*(NX + 1)*k + (NX + 1)*j + i;
281  elemToNode(ielem,1) = (NY + 1)*(NX + 1)*k + (NX + 1)*j + i + 1;
282  elemToNode(ielem,2) = (NY + 1)*(NX + 1)*k + (NX + 1)*(j + 1) + i + 1;
283  elemToNode(ielem,3) = (NY + 1)*(NX + 1)*k + (NX + 1)*(j + 1) + i;
284  elemToNode(ielem,4) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*j + i;
285  elemToNode(ielem,5) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*j + i + 1;
286  elemToNode(ielem,6) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*(j + 1) + i + 1;
287  elemToNode(ielem,7) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*(j + 1) + i;
288  ielem++;
289  }
290  }
291  }
292 
293  // Get edge connectivity
294  FieldContainer<int> edgeToNode(numEdges, numNodesPerEdge);
295  FieldContainer<int> elemToEdge(numElems, numEdgesPerElem);
296  int iedge = 0;
297  inode = 0;
298  for (int k=0; k<NZ+1; k++) {
299  for (int j=0; j<NY+1; j++) {
300  for (int i=0; i<NX+1; i++) {
301  if (i < NX){
302  edgeToNode(iedge,0) = inode;
303  edgeToNode(iedge,1) = inode + 1;
304  if (j < NY && k < NZ){
305  ielem=i+j*NX+k*NX*NY;
306  elemToEdge(ielem,0) = iedge;
307  if (j > 0)
308  elemToEdge(ielem-NX,2) = iedge;
309  if (k > 0)
310  elemToEdge(ielem-NX*NY,4) = iedge;
311  if (j > 0 && k > 0)
312  elemToEdge(ielem-NX*NY-NX,6) = iedge;
313  }
314  else if (j == NY && k == NZ){
315  ielem=i+(NY-1)*NX+(NZ-1)*NX*NY;
316  elemToEdge(ielem,6) = iedge;
317  }
318  else if (k == NZ && j < NY){
319  ielem=i+j*NX+(NZ-1)*NX*NY;
320  elemToEdge(ielem,4) = iedge;
321  if (j > 0)
322  elemToEdge(ielem-NX,6) = iedge;
323  }
324  else if (k < NZ && j == NY){
325  ielem=i+(NY-1)*NX+k*NX*NY;
326  elemToEdge(ielem,2) = iedge;
327  if (k > 0)
328  elemToEdge(ielem-NX*NY,6) = iedge;
329  }
330  iedge++;
331  }
332  if (j < NY){
333  edgeToNode(iedge,0) = inode;
334  edgeToNode(iedge,1) = inode + NX+1;
335  if (i < NX && k < NZ){
336  ielem=i+j*NX+k*NX*NY;
337  elemToEdge(ielem,3) = iedge;
338  if (i > 0)
339  elemToEdge(ielem-1,1) = iedge;
340  if (k > 0)
341  elemToEdge(ielem-NX*NY,7) = iedge;
342  if (i > 0 && k > 0)
343  elemToEdge(ielem-NX*NY-1,5) = iedge;
344  }
345  else if (i == NX && k == NZ){
346  ielem=NX-1+j*NX+(NZ-1)*NX*NY;
347  elemToEdge(ielem,5) = iedge;
348  }
349  else if (k == NZ && i < NX){
350  ielem=i+j*NX+(NZ-1)*NX*NY;
351  elemToEdge(ielem,7) = iedge;
352  if (i > 0)
353  elemToEdge(ielem-1,5) = iedge;
354  }
355  else if (k < NZ && i == NX){
356  ielem=NX-1+j*NX+k*NX*NY;
357  elemToEdge(ielem,1) = iedge;
358  if (k > 0)
359  elemToEdge(ielem-NX*NY,5) = iedge;
360  }
361  iedge++;
362  }
363  if (k < NZ){
364  edgeToNode(iedge,0) = inode;
365  edgeToNode(iedge,1) = inode + (NX+1)*(NY+1);
366  if (i < NX && j < NY){
367  ielem=i+j*NX+k*NX*NY;
368  elemToEdge(ielem,8) = iedge;
369  if (i > 0)
370  elemToEdge(ielem-1,9) = iedge;
371  if (j > 0)
372  elemToEdge(ielem-NX,11) = iedge;
373  if (i > 0 && j > 0)
374  elemToEdge(ielem-NX-1,10) = iedge;
375  }
376  else if (i == NX && j == NY){
377  ielem=NX-1+(NY-1)*NX+k*NX*NY;
378  elemToEdge(ielem,10) = iedge;
379  }
380  else if (j == NY && i < NX){
381  ielem=i+(NY-1)*NX+k*NX*NY;
382  elemToEdge(ielem,11) = iedge;
383  if (i > 0)
384  elemToEdge(ielem-1,10) = iedge;
385  }
386  else if (j < NY && i == NX){
387  ielem=NX-1+j*NX+k*NX*NY;
388  elemToEdge(ielem,9) = iedge;
389  if (j > 0)
390  elemToEdge(ielem-NX,10) = iedge;
391  }
392  iedge++;
393  }
394  inode++;
395  }
396  }
397  }
398 
399  // Find boundary edges
400  FieldContainer<int> edgeOnBoundary(numEdges);
401  for (int i=0; i<numEdges; i++){
402  if (nodeOnBoundary(edgeToNode(i,0)) && nodeOnBoundary(edgeToNode(i,1))){
403  edgeOnBoundary(i)=1;
404  }
405  }
406 
407  // Get face connectivity
408  FieldContainer<int> faceToNode(numFaces, numNodesPerFace);
409  FieldContainer<int> elemToFace(numElems, numFacesPerElem);
410  FieldContainer<int> faceToEdge(numFaces, numEdgesPerFace);
411  int iface = 0;
412  inode = 0;
413  for (int k=0; k<NZ+1; k++) {
414  for (int j=0; j<NY+1; j++) {
415  for (int i=0; i<NX+1; i++) {
416  if (i < NX && k < NZ) {
417  faceToNode(iface,0)=inode;
418  faceToNode(iface,1)=inode + 1;
419  faceToNode(iface,2)=inode + (NX+1)*(NY+1)+1;
420  faceToNode(iface,3)=inode + (NX+1)*(NY+1);
421  if (j < NY) {
422  ielem=i+j*NX+k*NX*NY;
423  faceToEdge(iface,0)=elemToEdge(ielem,0);
424  faceToEdge(iface,1)=elemToEdge(ielem,9);
425  faceToEdge(iface,2)=elemToEdge(ielem,4);
426  faceToEdge(iface,3)=elemToEdge(ielem,8);
427  elemToFace(ielem,0)=iface;
428  if (j > 0) {
429  elemToFace(ielem-NX,2)=iface;
430  }
431  }
432  else if (j == NY) {
433  ielem=i+(NY-1)*NX+k*NX*NY;
434  faceToEdge(iface,0)=elemToEdge(ielem,2);
435  faceToEdge(iface,1)=elemToEdge(ielem,10);
436  faceToEdge(iface,2)=elemToEdge(ielem,6);
437  faceToEdge(iface,3)=elemToEdge(ielem,11);
438  elemToFace(ielem,2)=iface;
439  }
440  iface++;
441  }
442  if (j < NY && k < NZ) {
443  faceToNode(iface,0)=inode;
444  faceToNode(iface,1)=inode + NX+1;
445  faceToNode(iface,2)=inode + (NX+1)*(NY+1) + NX+1;
446  faceToNode(iface,3)=inode + (NX+1)*(NY+1);
447  if (i < NX) {
448  ielem=i+j*NX+k*NX*NY;
449  faceToEdge(iface,0)=elemToEdge(ielem,3);
450  faceToEdge(iface,1)=elemToEdge(ielem,11);
451  faceToEdge(iface,2)=elemToEdge(ielem,7);
452  faceToEdge(iface,3)=elemToEdge(ielem,8);
453  elemToFace(ielem,3)=iface;
454  if (i > 0) {
455  elemToFace(ielem-1,1)=iface;
456  }
457  }
458  else if (i == NX) {
459  ielem=NX-1+j*NX+k*NX*NY;
460  faceToEdge(iface,0)=elemToEdge(ielem,1);
461  faceToEdge(iface,1)=elemToEdge(ielem,10);
462  faceToEdge(iface,2)=elemToEdge(ielem,5);
463  faceToEdge(iface,3)=elemToEdge(ielem,9);
464  elemToFace(ielem,1)=iface;
465  }
466  iface++;
467  }
468  if (i < NX && j < NY) {
469  faceToNode(iface,0)=inode;
470  faceToNode(iface,1)=inode + 1;
471  faceToNode(iface,2)=inode + NX+2;
472  faceToNode(iface,3)=inode + NX+1;
473  if (k < NZ) {
474  ielem=i+j*NX+k*NX*NY;
475  faceToEdge(iface,0)=elemToEdge(ielem,0);
476  faceToEdge(iface,1)=elemToEdge(ielem,1);
477  faceToEdge(iface,2)=elemToEdge(ielem,2);
478  faceToEdge(iface,3)=elemToEdge(ielem,3);
479  elemToFace(ielem,4)=iface;
480  if (k > 0) {
481  elemToFace(ielem-NX*NY,5)=iface;
482  }
483  }
484  else if (k == NZ) {
485  ielem=i+j*NX+(NZ-1)*NX*NY;
486  faceToEdge(iface,0)=elemToEdge(ielem,4);
487  faceToEdge(iface,1)=elemToEdge(ielem,5);
488  faceToEdge(iface,2)=elemToEdge(ielem,6);
489  faceToEdge(iface,3)=elemToEdge(ielem,7);
490  elemToFace(ielem,5)=iface;
491  }
492  iface++;
493  }
494  inode++;
495  }
496  }
497  }
498 
499  // Find boundary faces
500  FieldContainer<int> faceOnBoundary(numFaces);
501  for (int i=0; i<numFaces; i++){
502  if (nodeOnBoundary(faceToNode(i,0)) && nodeOnBoundary(faceToNode(i,1))
503  && nodeOnBoundary(faceToNode(i,2)) && nodeOnBoundary(faceToNode(i,3))){
504  faceOnBoundary(i)=1;
505  }
506  }
507 
508 
509 #define DUMP_DATA
510 #ifdef DUMP_DATA
511  // Output connectivity
512  ofstream fe2nout("elem2node.dat");
513  ofstream fe2eout("elem2edge.dat");
514  for (int k=0; k<NZ; k++) {
515  for (int j=0; j<NY; j++) {
516  for (int i=0; i<NX; i++) {
517  int ielem = i + j * NX + k * NX * NY;
518  for (int m=0; m<numNodesPerElem; m++){
519  fe2nout << elemToNode(ielem,m) <<" ";
520  }
521  fe2nout <<"\n";
522  for (int l=0; l<numEdgesPerElem; l++) {
523  fe2eout << elemToEdge(ielem,l) << " ";
524  }
525  fe2eout << "\n";
526  }
527  }
528  }
529  fe2nout.close();
530  fe2eout.close();
531 #endif
532 
533 #ifdef DUMP_DATA_EXTRA
534  ofstream fed2nout("edge2node.dat");
535  for (int i=0; i<numEdges; i++) {
536  fed2nout << edgeToNode(i,0) <<" ";
537  fed2nout << edgeToNode(i,1) <<"\n";
538  }
539  fed2nout.close();
540 
541  ofstream fBnodeout("nodeOnBndy.dat");
542  ofstream fBedgeout("edgeOnBndy.dat");
543  for (int i=0; i<numNodes; i++) {
544  fBnodeout << nodeOnBoundary(i) <<"\n";
545  }
546  for (int i=0; i<numEdges; i++) {
547  fBedgeout << edgeOnBoundary(i) <<"\n";
548  }
549  fBnodeout.close();
550  fBedgeout.close();
551 #endif
552 
553  // Set material properties using undeformed grid assuming each element has only one value of mu
554  FieldContainer<double> muVal(numElems);
555  for (int k=0; k<NZ; k++) {
556  for (int j=0; j<NY; j++) {
557  for (int i=0; i<NX; i++) {
558  int ielem = i + j * NX + k * NX * NY;
559  double midElemX = nodeCoord(elemToNode(ielem,0),0) + hx/2.0;
560  double midElemY = nodeCoord(elemToNode(ielem,0),1) + hy/2.0;
561  double midElemZ = nodeCoord(elemToNode(ielem,0),2) + hz/2.0;
562  if ( (midElemX > mu1LeftX) && (midElemY > mu1LeftY) && (midElemZ > mu1LeftZ) &&
563  (midElemX <= mu1RightX) && (midElemY <= mu1RightY) && (midElemZ <= mu1RightZ) ){
564  muVal(ielem) = mu1;
565  }
566  else {
567  muVal(ielem) = mu2;
568  }
569  }
570  }
571  }
572 
573  // Perturb mesh coordinates (only interior nodes)
574  if (randomMesh){
575  for (int k=1; k<NZ; k++) {
576  for (int j=1; j<NY; j++) {
577  for (int i=1; i<NX; i++) {
578  int inode = i + j * (NX + 1) + k * (NX + 1) * (NY + 1);
579  // random numbers between -1.0 and 1.0
580  double rx = 2.0 * (double)rand()/RAND_MAX - 1.0;
581  double ry = 2.0 * (double)rand()/RAND_MAX - 1.0;
582  double rz = 2.0 * (double)rand()/RAND_MAX - 1.0;
583  // limit variation to 1/4 edge length
584  nodeCoord(inode,0) = nodeCoord(inode,0) + 0.125 * hx * rx;
585  nodeCoord(inode,1) = nodeCoord(inode,1) + 0.125 * hy * ry;
586  nodeCoord(inode,2) = nodeCoord(inode,2) + 0.125 * hz * rz;
587  }
588  }
589  }
590  }
591 
592 #ifdef DUMP_DATA
593  // Print nodal coords
594  ofstream fcoordout("coords.dat");
595  for (int i=0; i<numNodes; i++) {
596  fcoordout << nodeCoord(i,0) <<" ";
597  fcoordout << nodeCoord(i,1) <<" ";
598  fcoordout << nodeCoord(i,2) <<"\n";
599  }
600  fcoordout.close();
601 #endif
602 
603 
604 // **************************** INCIDENCE MATRIX **************************************
605 
606  // Node to edge incidence matrix
607  *outStream << "Building incidence matrix ... \n\n";
608 
609  Epetra_SerialComm Comm;
610  Epetra_Map globalMapC(numEdges, 0, Comm);
611  Epetra_Map globalMapG(numNodes, 0, Comm);
612  Epetra_FECrsMatrix DGrad(Copy, globalMapC, globalMapG, 2);
613 
614  double vals[2];
615  vals[0]=-1.0; vals[1]=1.0;
616  for (int j=0; j<numEdges; j++){
617  int rowNum = j;
618  int colNum[2];
619  colNum[0] = edgeToNode(j,0);
620  colNum[1] = edgeToNode(j,1);
621  DGrad.InsertGlobalValues(1, &rowNum, 2, colNum, vals);
622  }
623 
624 
625 // ************************************ CUBATURE **************************************
626 
627  // Get numerical integration points and weights for element
628  *outStream << "Getting cubature ... \n\n";
629 
630  DefaultCubatureFactory<double> cubFactory;
631  int cubDegree = 2;
632  Teuchos::RCP<Cubature<double> > hexCub = cubFactory.create(hex_8, cubDegree);
633 
634  int cubDim = hexCub->getDimension();
635  int numCubPoints = hexCub->getNumPoints();
636 
637  FieldContainer<double> cubPoints(numCubPoints, cubDim);
638  FieldContainer<double> cubWeights(numCubPoints);
639 
640  hexCub->getCubature(cubPoints, cubWeights);
641 
642 
643  // Get numerical integration points and weights for hexahedron face
644  // (needed for rhs boundary term)
645 
646  // Define topology of the face parametrization domain as [-1,1]x[-1,1]
647  CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
648 
649  // Define cubature
650  DefaultCubatureFactory<double> cubFactoryFace;
651  Teuchos::RCP<Cubature<double> > hexFaceCubature = cubFactoryFace.create(paramQuadFace, 3);
652  int cubFaceDim = hexFaceCubature -> getDimension();
653  int numFacePoints = hexFaceCubature -> getNumPoints();
654 
655  // Define storage for cubature points and weights on [-1,1]x[-1,1]
656  FieldContainer<double> paramGaussWeights(numFacePoints);
657  FieldContainer<double> paramGaussPoints(numFacePoints,cubFaceDim);
658 
659  // Define storage for cubature points on workset faces
660  hexFaceCubature -> getCubature(paramGaussPoints, paramGaussWeights);
661 
662 
663 // ************************************** BASIS ***************************************
664 
665  // Define basis
666  *outStream << "Getting basis ... \n\n";
667  Basis_HCURL_HEX_I1_FEM<double, FieldContainer<double> > hexHCurlBasis;
668  Basis_HGRAD_HEX_C1_FEM<double, FieldContainer<double> > hexHGradBasis;
669 
670  int numFieldsC = hexHCurlBasis.getCardinality();
671  int numFieldsG = hexHGradBasis.getCardinality();
672 
673  // Evaluate basis at cubature points
674  FieldContainer<double> hexGVals(numFieldsG, numCubPoints);
675  FieldContainer<double> hexCVals(numFieldsC, numCubPoints, spaceDim);
676  FieldContainer<double> hexCurls(numFieldsC, numCubPoints, spaceDim);
677  FieldContainer<double> worksetCVals(numFieldsC, numFacePoints, spaceDim);
678 
679 
680  hexHCurlBasis.getValues(hexCVals, cubPoints, OPERATOR_VALUE);
681  hexHCurlBasis.getValues(hexCurls, cubPoints, OPERATOR_CURL);
682  hexHGradBasis.getValues(hexGVals, cubPoints, OPERATOR_VALUE);
683 
684 
685 // ******** LOOP OVER ELEMENTS TO CREATE LOCAL MASS and STIFFNESS MATRICES *************
686 
687 
688  *outStream << "Building mass and stiffness matrices ... \n\n";
689 
690  // Settings and data structures for mass and stiffness matrices
691  typedef CellTools<double> CellTools;
692  typedef FunctionSpaceTools fst;
693  //typedef ArrayTools art;
694  int numCells = 1;
695 
696  // Containers for nodes and edge signs
697  FieldContainer<double> hexNodes(numCells, numNodesPerElem, spaceDim);
698  FieldContainer<double> hexEdgeSigns(numCells, numFieldsC);
699  // Containers for Jacobian
700  FieldContainer<double> hexJacobian(numCells, numCubPoints, spaceDim, spaceDim);
701  FieldContainer<double> hexJacobInv(numCells, numCubPoints, spaceDim, spaceDim);
702  FieldContainer<double> hexJacobDet(numCells, numCubPoints);
703  // Containers for element HGRAD mass matrix
704  FieldContainer<double> massMatrixG(numCells, numFieldsG, numFieldsG);
705  FieldContainer<double> weightedMeasure(numCells, numCubPoints);
706  FieldContainer<double> weightedMeasureMuInv(numCells, numCubPoints);
707  FieldContainer<double> hexGValsTransformed(numCells, numFieldsG, numCubPoints);
708  FieldContainer<double> hexGValsTransformedWeighted(numCells, numFieldsG, numCubPoints);
709  // Containers for element HCURL mass matrix
710  FieldContainer<double> massMatrixC(numCells, numFieldsC, numFieldsC);
711  FieldContainer<double> hexCValsTransformed(numCells, numFieldsC, numCubPoints, spaceDim);
712  FieldContainer<double> hexCValsTransformedWeighted(numCells, numFieldsC, numCubPoints, spaceDim);
713  // Containers for element HCURL stiffness matrix
714  FieldContainer<double> stiffMatrixC(numCells, numFieldsC, numFieldsC);
715  FieldContainer<double> weightedMeasureMu(numCells, numCubPoints);
716  FieldContainer<double> hexCurlsTransformed(numCells, numFieldsC, numCubPoints, spaceDim);
717  FieldContainer<double> hexCurlsTransformedWeighted(numCells, numFieldsC, numCubPoints, spaceDim);
718  // Containers for right hand side vectors
719  FieldContainer<double> rhsDatag(numCells, numCubPoints, cubDim);
720  FieldContainer<double> rhsDatah(numCells, numCubPoints, cubDim);
721  FieldContainer<double> gC(numCells, numFieldsC);
722  FieldContainer<double> hC(numCells, numFieldsC);
723  FieldContainer<double> hCBoundary(numCells, numFieldsC);
724  FieldContainer<double> refGaussPoints(numFacePoints,spaceDim);
725  FieldContainer<double> worksetGaussPoints(numCells,numFacePoints,spaceDim);
726  FieldContainer<double> worksetJacobians(numCells, numFacePoints, spaceDim, spaceDim);
727  FieldContainer<double> worksetJacobInv(numCells, numFacePoints, spaceDim, spaceDim);
728  FieldContainer<double> worksetFaceN(numCells, numFacePoints, spaceDim);
729  FieldContainer<double> worksetFaceNweighted(numCells, numFacePoints, spaceDim);
730  FieldContainer<double> worksetVFieldVals(numCells, numFacePoints, spaceDim);
731  FieldContainer<double> worksetCValsTransformed(numCells, numFieldsC, numFacePoints, spaceDim);
732  FieldContainer<double> divuFace(numCells, numFacePoints);
733  FieldContainer<double> worksetFieldDotNormal(numCells, numFieldsC, numFacePoints);
734  // Container for cubature points in physical space
735  FieldContainer<double> physCubPoints(numCells,numCubPoints, cubDim);
736 
737 
738  // Global arrays in Epetra format
739  Epetra_FECrsMatrix MassG(Copy, globalMapG, numFieldsG);
740  Epetra_FECrsMatrix MassC(Copy, globalMapC, numFieldsC);
741  Epetra_FECrsMatrix StiffC(Copy, globalMapC, numFieldsC);
742  Epetra_FEVector rhsC(globalMapC);
743 
744 #ifdef DUMP_DATA
745  ofstream fSignsout("edgeSigns.dat");
746 #endif
747 
748  // *** Element loop ***
749  for (int k=0; k<numElems; k++) {
750 
751  // Physical cell coordinates
752  for (int i=0; i<numNodesPerElem; i++) {
753  hexNodes(0,i,0) = nodeCoord(elemToNode(k,i),0);
754  hexNodes(0,i,1) = nodeCoord(elemToNode(k,i),1);
755  hexNodes(0,i,2) = nodeCoord(elemToNode(k,i),2);
756  }
757 
758  // Edge signs
759  for (int j=0; j<numEdgesPerElem; j++) {
760  if (elemToNode(k,refEdgeToNode(j,0))==edgeToNode(elemToEdge(k,j),0) &&
761  elemToNode(k,refEdgeToNode(j,1))==edgeToNode(elemToEdge(k,j),1))
762  hexEdgeSigns(0,j) = 1.0;
763  else
764  hexEdgeSigns(0,j) = -1.0;
765 #ifdef DUMP_DATA
766  fSignsout << hexEdgeSigns(0,j) << " ";
767 #endif
768  }
769 #ifdef DUMP_DATA
770  fSignsout << "\n";
771 #endif
772 
773  // Compute cell Jacobians, their inverses and their determinants
774  CellTools::setJacobian(hexJacobian, cubPoints, hexNodes, hex_8);
775  CellTools::setJacobianInv(hexJacobInv, hexJacobian );
776  CellTools::setJacobianDet(hexJacobDet, hexJacobian );
777 
778 // ************************** Compute element HGrad mass matrices *******************************
779 
780  // transform to physical coordinates
781  fst::HGRADtransformVALUE<double>(hexGValsTransformed, hexGVals);
782 
783  // compute weighted measure
784  fst::computeCellMeasure<double>(weightedMeasure, hexJacobDet, cubWeights);
785 
786  // combine mu value with weighted measure
787  for (int nC = 0; nC < numCells; nC++){
788  for (int nPt = 0; nPt < numCubPoints; nPt++){
789  weightedMeasureMuInv(nC,nPt) = weightedMeasure(nC,nPt) / muVal(k);
790  }
791  }
792 
793  // multiply values with weighted measure
794  fst::multiplyMeasure<double>(hexGValsTransformedWeighted,
795  weightedMeasureMuInv, hexGValsTransformed);
796 
797  // integrate to compute element mass matrix
798  fst::integrate<double>(massMatrixG,
799  hexGValsTransformed, hexGValsTransformedWeighted, COMP_CPP);
800 
801  // assemble into global matrix
802  // int err = 0;
803  for (int row = 0; row < numFieldsG; row++){
804  for (int col = 0; col < numFieldsG; col++){
805  int rowIndex = elemToNode(k,row);
806  int colIndex = elemToNode(k,col);
807  double val = massMatrixG(0,row,col);
808  MassG.InsertGlobalValues(1, &rowIndex, 1, &colIndex, &val);
809  }
810  }
811 
812 // ************************** Compute element HCurl mass matrices *******************************
813 
814  // transform to physical coordinates
815  fst::HCURLtransformVALUE<double>(hexCValsTransformed, hexJacobInv,
816  hexCVals);
817 
818  // multiply by weighted measure
819  fst::multiplyMeasure<double>(hexCValsTransformedWeighted,
820  weightedMeasure, hexCValsTransformed);
821 
822  // integrate to compute element mass matrix
823  fst::integrate<double>(massMatrixC,
824  hexCValsTransformed, hexCValsTransformedWeighted,
825  COMP_CPP);
826 
827  // apply edge signs
828  fst::applyLeftFieldSigns<double>(massMatrixC, hexEdgeSigns);
829  fst::applyRightFieldSigns<double>(massMatrixC, hexEdgeSigns);
830 
831 
832  // assemble into global matrix
833  //err = 0;
834  for (int row = 0; row < numFieldsC; row++){
835  for (int col = 0; col < numFieldsC; col++){
836  int rowIndex = elemToEdge(k,row);
837  int colIndex = elemToEdge(k,col);
838  double val = massMatrixC(0,row,col);
839  MassC.InsertGlobalValues(1, &rowIndex, 1, &colIndex, &val);
840  }
841  }
842 
843 // ************************ Compute element HCurl stiffness matrices *****************************
844 
845  // transform to physical coordinates
846  fst::HCURLtransformCURL<double>(hexCurlsTransformed, hexJacobian, hexJacobDet,
847  hexCurls);
848 
849  // combine mu value with weighted measure
850  for (int nC = 0; nC < numCells; nC++){
851  for (int nPt = 0; nPt < numCubPoints; nPt++){
852  weightedMeasureMu(nC,nPt) = weightedMeasure(nC,nPt) / muVal(k);
853  }
854  }
855 
856  // multiply by weighted measure
857  fst::multiplyMeasure<double>(hexCurlsTransformedWeighted,
858  weightedMeasureMu, hexCurlsTransformed);
859 
860  // integrate to compute element stiffness matrix
861  fst::integrate<double>(stiffMatrixC,
862  hexCurlsTransformed, hexCurlsTransformedWeighted,
863  COMP_CPP);
864 
865  // apply edge signs
866  fst::applyLeftFieldSigns<double>(stiffMatrixC, hexEdgeSigns);
867  fst::applyRightFieldSigns<double>(stiffMatrixC, hexEdgeSigns);
868 
869  // assemble into global matrix
870  //err = 0;
871  for (int row = 0; row < numFieldsC; row++){
872  for (int col = 0; col < numFieldsC; col++){
873  int rowIndex = elemToEdge(k,row);
874  int colIndex = elemToEdge(k,col);
875  double val = stiffMatrixC(0,row,col);
876  StiffC.InsertGlobalValues(1, &rowIndex, 1, &colIndex, &val);
877  }
878  }
879 
880 // ******************************* Build right hand side ************************************
881 
882  // transform integration points to physical points
883  FieldContainer<double> physCubPoints(numCells,numCubPoints, cubDim);
884  CellTools::mapToPhysicalFrame(physCubPoints, cubPoints, hexNodes, hex_8);
885 
886  // evaluate right hand side functions at physical points
887  FieldContainer<double> rhsDatag(numCells, numCubPoints, cubDim);
888  FieldContainer<double> rhsDatah(numCells, numCubPoints, cubDim);
889  for (int nPt = 0; nPt < numCubPoints; nPt++){
890 
891  double x = physCubPoints(0,nPt,0);
892  double y = physCubPoints(0,nPt,1);
893  double z = physCubPoints(0,nPt,2);
894  double du1, du2, du3;
895 
896  evalCurlu(du1, du2, du3, x, y, z);
897  rhsDatag(0,nPt,0) = du1;
898  rhsDatag(0,nPt,1) = du2;
899  rhsDatag(0,nPt,2) = du3;
900 
901  evalGradDivu(du1, du2, du3, x, y, z);
902  rhsDatah(0,nPt,0) = du1;
903  rhsDatah(0,nPt,1) = du2;
904  rhsDatah(0,nPt,2) = du3;
905  }
906 
907  // integrate (g,curl w) term
908  fst::integrate<double>(gC, rhsDatag, hexCurlsTransformedWeighted,
909  COMP_CPP);
910 
911  // integrate -(grad h, w)
912  fst::integrate<double>(hC, rhsDatah, hexCValsTransformedWeighted,
913  COMP_CPP);
914 
915  // apply signs
916  fst::applyFieldSigns<double>(gC, hexEdgeSigns);
917  fst::applyFieldSigns<double>(hC, hexEdgeSigns);
918 
919 
920  // calculate boundary term, (h*w, n)_{\Gamma}
921  for (int i = 0; i < numFacesPerElem; i++){
922  if (faceOnBoundary(elemToFace(k,i))){
923 
924  // Map Gauss points on quad to reference face: paramGaussPoints -> refGaussPoints
925  CellTools::mapToReferenceSubcell(refGaussPoints,
926  paramGaussPoints,
927  2, i, hex_8);
928 
929  // Get basis values at points on reference cell
930  hexHCurlBasis.getValues(worksetCVals, refGaussPoints, OPERATOR_VALUE);
931 
932  // Compute Jacobians at Gauss pts. on reference face for all parent cells
933  CellTools::setJacobian(worksetJacobians,
934  refGaussPoints,
935  hexNodes, hex_8);
936  CellTools::setJacobianInv(worksetJacobInv, worksetJacobians );
937 
938  // transform to physical coordinates
939  fst::HCURLtransformVALUE<double>(worksetCValsTransformed, worksetJacobInv,
940  worksetCVals);
941 
942  // Map Gauss points on quad from ref. face to face workset: refGaussPoints -> worksetGaussPoints
943  CellTools::mapToPhysicalFrame(worksetGaussPoints,
944  refGaussPoints,
945  hexNodes, hex_8);
946 
947  // Compute face normals
948  CellTools::getPhysicalFaceNormals(worksetFaceN,
949  worksetJacobians,
950  i, hex_8);
951 
952  // multiply with weights
953  for(int nPt = 0; nPt < numFacePoints; nPt++){
954  for (int dim = 0; dim < spaceDim; dim++){
955  worksetFaceNweighted(0,nPt,dim) = worksetFaceN(0,nPt,dim) * paramGaussWeights(nPt);
956  } //dim
957  } //nPt
958 
959  fst::dotMultiplyDataField<double>(worksetFieldDotNormal, worksetFaceNweighted,
960  worksetCValsTransformed);
961 
962  // Evaluate div u at face points
963  for(int nPt = 0; nPt < numFacePoints; nPt++){
964 
965  double x = worksetGaussPoints(0, nPt, 0);
966  double y = worksetGaussPoints(0, nPt, 1);
967  double z = worksetGaussPoints(0, nPt, 2);
968 
969  divuFace(0,nPt)=evalDivu(x, y, z);
970  }
971 
972  // Integrate
973  fst::integrate<double>(hCBoundary, divuFace, worksetFieldDotNormal,
974  COMP_CPP);
975 
976  // apply signs
977  fst::applyFieldSigns<double>(hCBoundary, hexEdgeSigns);
978 
979  // add into hC term
980  for (int nF = 0; nF < numFieldsC; nF++){
981  hC(0,nF) = hC(0,nF) - hCBoundary(0,nF);
982  }
983 
984  } // if faceOnBoundary
985  } // numFaces
986 
987 
988  // assemble into global vector
989  for (int row = 0; row < numFieldsC; row++){
990  int rowIndex = elemToEdge(k,row);
991  double val = gC(0,row)-hC(0,row);
992  rhsC.SumIntoGlobalValues(1, &rowIndex, &val);
993  }
994 
995 
996  } // *** end element loop ***
997 
998 #ifdef DUMP_DATA
999  fSignsout.close();
1000 #endif
1001 
1002  // Assemble over multiple processors, if necessary
1003  MassG.GlobalAssemble(); MassG.FillComplete();
1004  MassC.GlobalAssemble(); MassC.FillComplete();
1005  StiffC.GlobalAssemble(); StiffC.FillComplete();
1006  rhsC.GlobalAssemble();
1007  DGrad.GlobalAssemble(); DGrad.FillComplete(MassG.RowMap(),MassC.RowMap());
1008 
1009 
1010  // Build the inverse diagonal for MassG
1011  Epetra_CrsMatrix MassGinv(Copy,MassG.RowMap(),MassG.RowMap(),1);
1012  Epetra_Vector DiagG(MassG.RowMap());
1013 
1014  DiagG.PutScalar(1.0);
1015  MassG.Multiply(false,DiagG,DiagG);
1016  for(int i=0; i<DiagG.MyLength(); i++) {
1017  DiagG[i]=1.0/DiagG[i];
1018  }
1019  for(int i=0; i<DiagG.MyLength(); i++) {
1020  int CID=MassG.GCID(i);
1021  MassGinv.InsertGlobalValues(MassG.GRID(i),1,&(DiagG[i]),&CID);
1022  }
1023  MassGinv.FillComplete();
1024 
1025  // Set value to zero on diagonal that cooresponds to boundary node
1026  for(int i=0;i<numNodes;i++) {
1027  if (nodeOnBoundary(i)){
1028  double val=0.0;
1029  MassGinv.ReplaceGlobalValues(i,1,&val,&i);
1030  }
1031  }
1032 
1033  int numEntries;
1034  double *values;
1035  int *cols;
1036 
1037  // Adjust matrices and rhs due to boundary conditions
1038  for (int row = 0; row<numEdges; row++){
1039  MassC.ExtractMyRowView(row,numEntries,values,cols);
1040  for (int i=0; i<numEntries; i++){
1041  if (edgeOnBoundary(cols[i])) {
1042  values[i]=0;
1043  }
1044  }
1045  StiffC.ExtractMyRowView(row,numEntries,values,cols);
1046  for (int i=0; i<numEntries; i++){
1047  if (edgeOnBoundary(cols[i])) {
1048  values[i]=0;
1049  }
1050  }
1051  }
1052  for (int row = 0; row<numEdges; row++){
1053  if (edgeOnBoundary(row)) {
1054  int rowindex = row;
1055  StiffC.ExtractMyRowView(row,numEntries,values,cols);
1056  for (int i=0; i<numEntries; i++){
1057  values[i]=0;
1058  }
1059  MassC.ExtractMyRowView(row,numEntries,values,cols);
1060  for (int i=0; i<numEntries; i++){
1061  values[i]=0;
1062  }
1063  rhsC[0][row]=0;
1064  double val = 1.0;
1065  StiffC.ReplaceGlobalValues(1, &rowindex, 1, &rowindex, &val);
1066  }
1067  }
1068 
1069 
1070 #ifdef DUMP_DATA
1071  // Dump matrices to disk
1072  EpetraExt::RowMatrixToMatlabFile("mag_m0inv_matrix.dat",MassGinv);
1073  EpetraExt::RowMatrixToMatlabFile("mag_m1_matrix.dat",MassC);
1074  EpetraExt::RowMatrixToMatlabFile("mag_k1_matrix.dat",StiffC);
1075  EpetraExt::RowMatrixToMatlabFile("mag_t0_matrix.dat",DGrad);
1076  EpetraExt::MultiVectorToMatrixMarketFile("mag_rhs1_vector.dat",rhsC,0,0,false);
1077  fSignsout.close();
1078 #endif
1079 
1080 
1081  std::cout << "End Result: TEST PASSED\n";
1082 
1083  // reset format state of std::cout
1084  std::cout.copyfmt(oldFormatState);
1085 
1086  return 0;
1087 }
1088 
1089 // Calculates value of exact solution u
1090  int evalu(double & uExact0, double & uExact1, double & uExact2, double & x, double & y, double & z)
1091  {
1092  // function 1
1093  uExact0 = exp(x+y+z)*(y+1.0)*(y-1.0)*(z+1.0)*(z-1.0);
1094  uExact1 = exp(x+y+z)*(x+1.0)*(x-1.0)*(z+1.0)*(z-1.0);
1095  uExact2 = exp(x+y+z)*(x+1.0)*(x-1.0)*(y+1.0)*(y-1.0);
1096 
1097  /*
1098  // function 2
1099  uExact0 = cos(M_PI*x)*exp(y*z)*(y+1.0)*(y-1.0)*(z+1.0)*(z-1.0);
1100  uExact1 = cos(M_PI*y)*exp(x*z)*(x+1.0)*(x-1.0)*(z+1.0)*(z-1.0);
1101  uExact2 = cos(M_PI*z)*exp(x*y)*(x+1.0)*(x-1.0)*(y+1.0)*(y-1.0);
1102 
1103  // function 3
1104  uExact0 = cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
1105  uExact1 = sin(M_PI*x)*cos(M_PI*y)*sin(M_PI*z);
1106  uExact2 = sin(M_PI*x)*sin(M_PI*y)*cos(M_PI*z);
1107 
1108  // function 4
1109  uExact0 = (y*y - 1.0)*(z*z-1.0);
1110  uExact1 = (x*x - 1.0)*(z*z-1.0);
1111  uExact2 = (x*x - 1.0)*(y*y-1.0);
1112  */
1113 
1114  return 0;
1115  }
1116 
1117 // Calculates divergence of exact solution u
1118  double evalDivu(double & x, double & y, double & z)
1119  {
1120  // function 1
1121  double divu = exp(x+y+z)*(y+1.0)*(y-1.0)*(z+1.0)*(z-1.0)
1122  + exp(x+y+z)*(x+1.0)*(x-1.0)*(z+1.0)*(z-1.0)
1123  + exp(x+y+z)*(x+1.0)*(x-1.0)*(y+1.0)*(y-1.0);
1124 
1125  /*
1126  // function 2
1127  double divu = -M_PI*sin(M_PI*x)*exp(y*z)*(y+1.0)*(y-1.0)*(z+1.0)*(z-1.0)
1128  -M_PI*sin(M_PI*y)*exp(x*z)*(x+1.0)*(x-1.0)*(z+1.0)*(z-1.0)
1129  -M_PI*sin(M_PI*z)*exp(x*y)*(x+1.0)*(x-1.0)*(y+1.0)*(y-1.0);
1130 
1131  // function 3
1132  double divu = -3.0*M_PI*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
1133 
1134  // function 4
1135  double divu = 0.0;
1136  */
1137 
1138  return divu;
1139  }
1140 
1141 
1142 // Calculates curl of exact solution u
1143  int evalCurlu(double & curlu0, double & curlu1, double & curlu2, double & x, double & y, double & z)
1144  {
1145 
1146  // function 1
1147  double duxdy = exp(x+y+z)*(z*z-1.0)*(y*y+2.0*y-1.0);
1148  double duxdz = exp(x+y+z)*(y*y-1.0)*(z*z+2.0*z-1.0);
1149  double duydx = exp(x+y+z)*(z*z-1.0)*(x*x+2.0*x-1.0);
1150  double duydz = exp(x+y+z)*(x*x-1.0)*(z*z+2.0*z-1.0);
1151  double duzdx = exp(x+y+z)*(y*y-1.0)*(x*x+2.0*x-1.0);
1152  double duzdy = exp(x+y+z)*(x*x-1.0)*(y*y+2.0*y-1.0);
1153 
1154  /*
1155  // function 2
1156  double duxdy = cos(M_PI*x)*exp(y*z)*(z+1.0)*(z-1.0)*(z*(y+1.0)*(y-1.0) + 2.0*y);
1157  double duxdz = cos(M_PI*x)*exp(y*z)*(y+1.0)*(y-1.0)*(y*(z+1.0)*(z-1.0) + 2.0*z);
1158  double duydx = cos(M_PI*y)*exp(x*z)*(z+1.0)*(z-1.0)*(z*(x+1.0)*(x-1.0) + 2.0*x);
1159  double duydz = cos(M_PI*y)*exp(x*z)*(x+1.0)*(x-1.0)*(x*(z+1.0)*(z-1.0) + 2.0*z);
1160  double duzdx = cos(M_PI*z)*exp(x*y)*(y+1.0)*(y-1.0)*(y*(x+1.0)*(x-1.0) + 2.0*x);
1161  double duzdy = cos(M_PI*z)*exp(x*y)*(x+1.0)*(x-1.0)*(x*(y+1.0)*(y-1.0) + 2.0*y);
1162 
1163  // function 3
1164  double duxdy = M_PI*cos(M_PI*x)*cos(M_PI*y)*sin(M_PI*z);
1165  double duxdz = M_PI*cos(M_PI*x)*sin(M_PI*y)*cos(M_PI*z);
1166  double duydx = M_PI*cos(M_PI*x)*cos(M_PI*y)*sin(M_PI*z);
1167  double duydz = M_PI*sin(M_PI*x)*cos(M_PI*y)*cos(M_PI*z);
1168  double duzdx = M_PI*cos(M_PI*x)*sin(M_PI*y)*cos(M_PI*z);
1169  double duzdy = M_PI*sin(M_PI*x)*cos(M_PI*y)*cos(M_PI*z);
1170 
1171  // function 4
1172  double duxdy = 2.0*y*(z*z-1);
1173  double duxdz = 2.0*z*(y*y-1);
1174  double duydx = 2.0*x*(z*z-1);
1175  double duydz = 2.0*z*(x*x-1);
1176  double duzdx = 2.0*x*(y*y-1);
1177  double duzdy = 2.0*y*(x*x-1);
1178 */
1179 
1180  curlu0 = duzdy - duydz;
1181  curlu1 = duxdz - duzdx;
1182  curlu2 = duydx - duxdy;
1183 
1184  return 0;
1185  }
1186 
1187 // Calculates gradient of the divergence of exact solution u
1188  int evalGradDivu(double & gradDivu0, double & gradDivu1, double & gradDivu2, double & x, double & y, double & z)
1189 {
1190 
1191  // function 1
1192  gradDivu0 = exp(x+y+z)*((y*y-1.0)*(z*z-1.0)+(x*x+2.0*x-1.0)*(z*z-1.0)+(x*x+2.0*x-1.0)*(y*y-1.0));
1193  gradDivu1 = exp(x+y+z)*((y*y+2.0*y-1.0)*(z*z-1.0)+(x*x-1.0)*(z*z-1.0)+(x*x-1.0)*(y*y+2.0*y-1.0));
1194  gradDivu2 = exp(x+y+z)*((y*y-1.0)*(z*z+2.0*z-1.0)+(x*x-1.0)*(z*z+2.0*z-1.0)+(x*x-1.0)*(y*y-1.0));
1195 
1196  /*
1197  // function 2
1198  gradDivu0 = -M_PI*M_PI*cos(M_PI*x)*exp(y*z)*(y+1.0)*(y-1.0)*(z+1.0)*(z-1.0)
1199  -M_PI*sin(M_PI*y)*exp(x*z)*(z+1.0)*(z-1.0)*(z*(x+1.0)*(x-1.0)+2.0*x)
1200  -M_PI*sin(M_PI*z)*exp(x*y)*(y+1.0)*(y-1.0)*(y*(x+1.0)*(x-1.0)+2.0*x);
1201  gradDivu1 = -M_PI*sin(M_PI*x)*exp(y*z)*(z+1.0)*(z-1.0)*(z*(y+1.0)*(y-1.0)+2.0*y)
1202  -M_PI*M_PI*cos(M_PI*y)*exp(x*z)*(x+1.0)*(x-1.0)*(z+1.0)*(z-1.0)
1203  -M_PI*sin(M_PI*z)*exp(x*y)*(x+1.0)*(x-1.0)*(x*(y+1.0)*(y-1.0)+2.0*y);
1204  gradDivu2 = -M_PI*sin(M_PI*x)*exp(y*z)*(y+1.0)*(y-1.0)*(y*(z+1.0)*(z-1.0)+2.0*z)
1205  -M_PI*sin(M_PI*y)*exp(x*z)*(x+1.0)*(x-1.0)*(x*(z+1.0)*(z-1.0)+2.0*z)
1206  -M_PI*M_PI*cos(M_PI*z)*exp(x*y)*(x+1.0)*(x-1.0)*(y+1.0)*(y-1.0);
1207 
1208  // function 3
1209  gradDivu0 = -3.0*M_PI*M_PI*cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
1210  gradDivu1 = -3.0*M_PI*M_PI*sin(M_PI*x)*cos(M_PI*y)*sin(M_PI*z);
1211  gradDivu2 = -3.0*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y)*cos(M_PI*z);
1212 
1213  // function 4
1214  gradDivu0 = 0;
1215  gradDivu1 = 0;
1216  gradDivu2 = 0;
1217  */
1218 
1219  return 0;
1220 }
Intrepid::Basis::getCardinality
virtual int getCardinality() const
Returns cardinality of the basis.
Definition: Intrepid_BasisDef.hpp:100
Intrepid_CellTools.hpp
Header file for the Intrepid::CellTools class.
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid_ArrayTools.hpp
Header file for utility class to provide array tools, such as tensor contractions, etc.
Intrepid::FunctionSpaceTools
Defines expert-level interfaces for the evaluation of functions and operators in physical space (supp...
Definition: Intrepid_FunctionSpaceTools.hpp:69
Intrepid::Basis_HCURL_HEX_I1_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
Definition: Intrepid_HCURL_HEX_I1_FEMDef.hpp:100
Intrepid_Utils.hpp
Intrepid utilities.
Intrepid_DefaultCubatureFactory.hpp
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Intrepid::Basis_HCURL_HEX_I1_FEM
Implementation of the default H(curl)-compatible FEM basis of degree 1 on Hexahedron cell...
Definition: Intrepid_HCURL_HEX_I1_FEM.hpp:115
Intrepid::FieldContainer
Implementation of a templated lexicographical container for a multi-indexed scalar quantity...
Definition: Intrepid_FieldContainer.hpp:78
Intrepid_FunctionSpaceTools.hpp
Header file for the Intrepid::FunctionSpaceTools class.
Intrepid_RealSpaceTools.hpp
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
Intrepid::Basis_HGRAD_HEX_C1_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Hexahedron cell...
Definition: Intrepid_HGRAD_HEX_C1_FEM.hpp:91
Intrepid::DefaultCubatureFactory
A factory class that generates specific instances of cubatures.
Definition: Intrepid_DefaultCubatureFactory.hpp:77
Intrepid::DefaultCubatureFactory::create
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > &degree)
Factory method.
Definition: Intrepid_DefaultCubatureFactoryDef.hpp:53
Intrepid::Basis_HGRAD_HEX_C1_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
Definition: Intrepid_HGRAD_HEX_C1_FEMDef.hpp:100
Intrepid::CellTools
A stateless class for operations on cell data. Provides methods for:
Definition: Intrepid_CellTools.hpp:109
Intrepid_HCURL_HEX_I1_FEM.hpp
Header file for the Intrepid::HCURL_HEX_I1_FEM class.
Generated by   doxygen 1.8.5