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