Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/example/Drivers/example_01.cpp
00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //                           Intrepid Package
00005 //                 Copyright (2007) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00092 // Intrepid includes
00093 #include "Intrepid_FunctionSpaceTools.hpp"
00094 #include "Intrepid_FieldContainer.hpp"
00095 #include "Intrepid_CellTools.hpp"
00096 #include "Intrepid_ArrayTools.hpp"
00097 #include "Intrepid_HCURL_HEX_I1_FEM.hpp"
00098 #include "Intrepid_HGRAD_HEX_C1_FEM.hpp"
00099 #include "Intrepid_RealSpaceTools.hpp"
00100 #include "Intrepid_DefaultCubatureFactory.hpp"
00101 #include "Intrepid_Utils.hpp"
00102 
00103 // Epetra includes
00104 #include "Epetra_Time.h"
00105 #include "Epetra_Map.h"
00106 #include "Epetra_SerialComm.h"
00107 #include "Epetra_FECrsMatrix.h"
00108 #include "Epetra_FEVector.h"
00109 #include "Epetra_Vector.h"
00110 
00111 // Teuchos includes
00112 #include "Teuchos_oblackholestream.hpp"
00113 #include "Teuchos_RCP.hpp"
00114 #include "Teuchos_BLAS.hpp"
00115 
00116 // Shards includes
00117 #include "Shards_CellTopology.hpp"
00118 
00119 // EpetraExt includes
00120 #include "EpetraExt_RowMatrixOut.h"
00121 #include "EpetraExt_MultiVectorOut.h"
00122 
00123 using namespace std;
00124 using namespace Intrepid;
00125 
00126 // Functions to evaluate exact solution and derivatives
00127 int evalu(double & uExact0, double & uExact1, double & uExact2, double & x, double & y, double & z);
00128 double evalDivu(double & x, double & y, double & z);
00129 int evalCurlu(double & curlu0, double & curlu1, double & curlu2, double & x, double & y, double & z);
00130 int evalGradDivu(double & gradDivu0, double & gradDivu1, double & gradDivu2, double & x, double & y, double & z);
00131 
00132 int main(int argc, char *argv[]) {
00133 
00134    //Check number of arguments
00135    if (argc < 13) {
00136       std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
00137       std::cout <<"Usage:\n\n";
00138       std::cout <<"  ./Intrepid_example_Drivers_Example_01.exe NX NY NZ randomMesh mu1 mu2 mu1LX mu1RX mu1LY mu1RY mu1LZ mu1RZ verbose\n\n";
00139       std::cout <<" where \n";
00140       std::cout <<"   int NX              - num intervals in x direction (assumed box domain, -1,1) \n";
00141       std::cout <<"   int NY              - num intervals in y direction (assumed box domain, -1,1) \n";
00142       std::cout <<"   int NZ              - num intervals in z direction (assumed box domain, -1,1) \n";
00143       std::cout <<"   int randomMesh      - 1 if mesh randomizer is to be used 0 if not \n";
00144       std::cout <<"   double mu1          - material property value for region 1 \n";
00145       std::cout <<"   double mu2          - material property value for region 2 \n";
00146       std::cout <<"   double mu1LX        - left X boundary for region 1 \n";
00147       std::cout <<"   double mu1RX        - right X boundary for region 1 \n";
00148       std::cout <<"   double mu1LY        - left Y boundary for region 1 \n";
00149       std::cout <<"   double mu1RY        - right Y boundary for region 1 \n";
00150       std::cout <<"   double mu1LZ        - bottom Z boundary for region 1 \n";
00151       std::cout <<"   double mu1RZ        - top Z boundary for region 1 \n";
00152       std::cout <<"   verbose (optional)  - any character, indicates verbose output \n\n";
00153       exit(1);
00154    }
00155   
00156   // This little trick lets us print to std::cout only if
00157   // a (dummy) command-line argument is provided.
00158   int iprint     = argc - 1;
00159   Teuchos::RCP<std::ostream> outStream;
00160   Teuchos::oblackholestream bhs; // outputs nothing
00161   if (iprint > 12)
00162     outStream = Teuchos::rcp(&std::cout, false);
00163   else
00164     outStream = Teuchos::rcp(&bhs, false);
00165   
00166   // Save the format state of the original std::cout.
00167   Teuchos::oblackholestream oldFormatState;
00168   oldFormatState.copyfmt(std::cout);
00169   
00170   *outStream \
00171     << "===============================================================================\n" \
00172     << "|                                                                             |\n" \
00173     << "|  Example: Generate Mass and Stiffness Matrices and Right-Hand Side Vector   |\n"
00174     << "|    for Div-Curl System on Hexahedral Mesh with Curl-Conforming Elements     |\n" \
00175     << "|                                                                             |\n" \
00176     << "|  Questions? Contact  Pavel Bochev  (pbboche@sandia.gov),                    |\n" \
00177     << "|                      Denis Ridzal  (dridzal@sandia.gov),                    |\n" \
00178     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00179     << "|                                                                             |\n" \
00180     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00181     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00182     << "|                                                                             |\n" \
00183     << "===============================================================================\n";
00184 
00185 
00186 // ************************************ GET INPUTS **************************************
00187 
00188   /* In the implementation for discontinuous material properties only the boundaries for
00189      region 1, associated with mu1, are input. The remainder of the grid is assumed to use mu2.
00190      Note that the material properties are assigned using the undeformed grid. */
00191 
00192     int NX            = atoi(argv[1]);  // num intervals in x direction (assumed box domain, -1,1)
00193     int NY            = atoi(argv[2]);  // num intervals in y direction (assumed box domain, -1,1)
00194     int NZ            = atoi(argv[3]);  // num intervals in z direction (assumed box domain, -1,1)
00195     int randomMesh    = atoi(argv[4]);  // 1 if mesh randomizer is to be used 0 if not
00196     double mu1        = atof(argv[5]);  // material property value for region 1
00197     double mu2        = atof(argv[6]);  // material property value for region 2
00198     double mu1LeftX   = atof(argv[7]);  // left X boundary for region 1
00199     double mu1RightX  = atof(argv[8]);  // right X boundary for region 1
00200     double mu1LeftY   = atof(argv[9]);  // left Y boundary for region 1
00201     double mu1RightY  = atof(argv[10]); // right Y boundary for region 1
00202     double mu1LeftZ   = atof(argv[11]); // left Z boundary for region 1
00203     double mu1RightZ  = atof(argv[12]); // right Z boundary for region 1
00204 
00205 // *********************************** CELL TOPOLOGY **********************************
00206 
00207    // Get cell topology for base hexahedron
00208     typedef shards::CellTopology    CellTopology;
00209     CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() );
00210 
00211    // Get dimensions 
00212     int numNodesPerElem = hex_8.getNodeCount();
00213     int numEdgesPerElem = hex_8.getEdgeCount();
00214     int numFacesPerElem = hex_8.getSideCount();
00215     int numNodesPerEdge = 2;
00216     int numNodesPerFace = 4;
00217     int numEdgesPerFace = 4;
00218     int spaceDim = hex_8.getDimension();
00219 
00220    // Build reference element edge to node map
00221     FieldContainer<int> refEdgeToNode(numEdgesPerElem,numNodesPerEdge);
00222     for (int i=0; i<numEdgesPerElem; i++){
00223         refEdgeToNode(i,0)=hex_8.getNodeMap(1, i, 0);
00224         refEdgeToNode(i,1)=hex_8.getNodeMap(1, i, 1);
00225     }
00226 
00227 // *********************************** GENERATE MESH ************************************
00228 
00229     *outStream << "Generating mesh ... \n\n";
00230 
00231     *outStream << "    NX" << "   NY" << "   NZ\n";
00232     *outStream << std::setw(5) << NX <<
00233                  std::setw(5) << NY <<
00234                  std::setw(5) << NZ << "\n\n";
00235 
00236    // Print mesh information  
00237     int numElems = NX*NY*NZ;
00238     int numNodes = (NX+1)*(NY+1)*(NZ+1);
00239     int numEdges = (NX)*(NY + 1)*(NZ + 1) + (NX + 1)*(NY)*(NZ + 1) + (NX + 1)*(NY + 1)*(NZ);
00240     int numFaces = (NX)*(NY)*(NZ + 1) + (NX)*(NY + 1)*(NZ) + (NX + 1)*(NY)*(NZ);
00241     *outStream << " Number of Elements: " << numElems << " \n";
00242     *outStream << "    Number of Nodes: " << numNodes << " \n";
00243     *outStream << "    Number of Edges: " << numEdges << " \n";
00244     *outStream << "    Number of Faces: " << numFaces << " \n\n";
00245 
00246    // Cube
00247     double leftX = -1.0, rightX = 1.0;
00248     double leftY = -1.0, rightY = 1.0;
00249     double leftZ = -1.0, rightZ = 1.0;
00250 
00251    // Mesh spacing
00252     double hx = (rightX-leftX)/((double)NX);
00253     double hy = (rightY-leftY)/((double)NY);
00254     double hz = (rightZ-leftZ)/((double)NZ);
00255 
00256     // Get nodal coordinates
00257     FieldContainer<double> nodeCoord(numNodes, spaceDim);
00258     FieldContainer<int> nodeOnBoundary(numNodes);
00259     int inode = 0;
00260     for (int k=0; k<NZ+1; k++) {
00261       for (int j=0; j<NY+1; j++) {
00262         for (int i=0; i<NX+1; i++) {
00263           nodeCoord(inode,0) = leftX + (double)i*hx;
00264           nodeCoord(inode,1) = leftY + (double)j*hy;
00265           nodeCoord(inode,2) = leftZ + (double)k*hz;
00266           if (k==0 || j==0 || i==0 || k==NZ || j==NY || i==NX){
00267              nodeOnBoundary(inode)=1; 
00268           }
00269           inode++;
00270         }
00271       }
00272     }
00273     
00274    // Element to Node map
00275     FieldContainer<int> elemToNode(numElems, numNodesPerElem);
00276     int ielem = 0;
00277     for (int k=0; k<NZ; k++) {
00278       for (int j=0; j<NY; j++) {
00279         for (int i=0; i<NX; i++) {
00280           elemToNode(ielem,0) = (NY + 1)*(NX + 1)*k + (NX + 1)*j + i;
00281           elemToNode(ielem,1) = (NY + 1)*(NX + 1)*k + (NX + 1)*j + i + 1;
00282           elemToNode(ielem,2) = (NY + 1)*(NX + 1)*k + (NX + 1)*(j + 1) + i + 1;
00283           elemToNode(ielem,3) = (NY + 1)*(NX + 1)*k + (NX + 1)*(j + 1) + i;
00284           elemToNode(ielem,4) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*j + i;
00285           elemToNode(ielem,5) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*j + i + 1;
00286           elemToNode(ielem,6) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*(j + 1) + i + 1;
00287           elemToNode(ielem,7) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*(j + 1) + i;
00288           ielem++;
00289         }
00290       }
00291     }
00292 
00293   // Get edge connectivity
00294     FieldContainer<int> edgeToNode(numEdges, numNodesPerEdge);
00295     FieldContainer<int> elemToEdge(numElems, numEdgesPerElem);
00296     int iedge = 0;
00297     inode = 0;
00298     for (int k=0; k<NZ+1; k++) {
00299       for (int j=0; j<NY+1; j++) {
00300         for (int i=0; i<NX+1; i++) {
00301            if (i < NX){
00302                edgeToNode(iedge,0) = inode;
00303                edgeToNode(iedge,1) = inode + 1;
00304                if (j < NY && k < NZ){
00305                   ielem=i+j*NX+k*NX*NY;
00306                   elemToEdge(ielem,0) = iedge;
00307                   if (j > 0)
00308                      elemToEdge(ielem-NX,2) = iedge; 
00309                   if (k > 0)
00310                      elemToEdge(ielem-NX*NY,4) = iedge; 
00311                   if (j > 0 && k > 0)
00312                      elemToEdge(ielem-NX*NY-NX,6) = iedge; 
00313                 }
00314                else if (j == NY && k == NZ){
00315                   ielem=i+(NY-1)*NX+(NZ-1)*NX*NY;
00316                   elemToEdge(ielem,6) = iedge;
00317                 }
00318                else if (k == NZ && j < NY){
00319                   ielem=i+j*NX+(NZ-1)*NX*NY;
00320                   elemToEdge(ielem,4) = iedge;
00321                   if (j > 0)
00322                     elemToEdge(ielem-NX,6) = iedge;
00323                 }
00324                else if (k < NZ && j == NY){
00325                   ielem=i+(NY-1)*NX+k*NX*NY;
00326                   elemToEdge(ielem,2) = iedge;
00327                   if (k > 0)
00328                      elemToEdge(ielem-NX*NY,6) = iedge;
00329                 }
00330                iedge++;
00331             }
00332            if (j < NY){
00333                edgeToNode(iedge,0) = inode;
00334                edgeToNode(iedge,1) = inode + NX+1;
00335                if (i < NX && k < NZ){
00336                   ielem=i+j*NX+k*NX*NY;
00337                   elemToEdge(ielem,3) = iedge;
00338                   if (i > 0)
00339                      elemToEdge(ielem-1,1) = iedge; 
00340                   if (k > 0)
00341                      elemToEdge(ielem-NX*NY,7) = iedge; 
00342                   if (i > 0 && k > 0)
00343                      elemToEdge(ielem-NX*NY-1,5) = iedge; 
00344                 }
00345                else if (i == NX && k == NZ){
00346                   ielem=NX-1+j*NX+(NZ-1)*NX*NY;
00347                   elemToEdge(ielem,5) = iedge;
00348                 }
00349                else if (k == NZ && i < NX){
00350                   ielem=i+j*NX+(NZ-1)*NX*NY;
00351                   elemToEdge(ielem,7) = iedge;
00352                   if (i > 0)
00353                     elemToEdge(ielem-1,5) = iedge;
00354                 }
00355                else if (k < NZ && i == NX){
00356                   ielem=NX-1+j*NX+k*NX*NY;
00357                   elemToEdge(ielem,1) = iedge;
00358                   if (k > 0)
00359                      elemToEdge(ielem-NX*NY,5) = iedge;
00360                 }
00361                iedge++;
00362             }
00363            if (k < NZ){
00364                edgeToNode(iedge,0) = inode;
00365                edgeToNode(iedge,1) = inode + (NX+1)*(NY+1);
00366                if (i < NX && j < NY){
00367                   ielem=i+j*NX+k*NX*NY;
00368                   elemToEdge(ielem,8) = iedge;
00369                   if (i > 0)
00370                      elemToEdge(ielem-1,9) = iedge; 
00371                   if (j > 0)
00372                      elemToEdge(ielem-NX,11) = iedge; 
00373                   if (i > 0 && j > 0)
00374                      elemToEdge(ielem-NX-1,10) = iedge; 
00375                 }
00376                else if (i == NX && j == NY){
00377                   ielem=NX-1+(NY-1)*NX+k*NX*NY;
00378                   elemToEdge(ielem,10) = iedge;
00379                 }
00380                else if (j == NY && i < NX){
00381                   ielem=i+(NY-1)*NX+k*NX*NY;
00382                   elemToEdge(ielem,11) = iedge;
00383                   if (i > 0)
00384                     elemToEdge(ielem-1,10) = iedge;
00385                 }
00386                else if (j < NY && i == NX){
00387                   ielem=NX-1+j*NX+k*NX*NY;
00388                   elemToEdge(ielem,9) = iedge;
00389                   if (j > 0)
00390                      elemToEdge(ielem-NX,10) = iedge;
00391                 }
00392                iedge++;
00393             }
00394             inode++;
00395          }
00396       }
00397    }
00398 
00399    // Find boundary edges  
00400     FieldContainer<int> edgeOnBoundary(numEdges);
00401     for (int i=0; i<numEdges; i++){
00402        if (nodeOnBoundary(edgeToNode(i,0)) && nodeOnBoundary(edgeToNode(i,1))){
00403            edgeOnBoundary(i)=1;
00404        }
00405     }
00406 
00407    // Get face connectivity
00408     FieldContainer<int> faceToNode(numFaces, numNodesPerFace);
00409     FieldContainer<int> elemToFace(numElems, numFacesPerElem);
00410     FieldContainer<int> faceToEdge(numFaces, numEdgesPerFace);
00411     int iface = 0;
00412     inode = 0;
00413     for (int k=0; k<NZ+1; k++) {
00414       for (int j=0; j<NY+1; j++) {
00415         for (int i=0; i<NX+1; i++) {
00416            if (i < NX && k < NZ) {
00417               faceToNode(iface,0)=inode;
00418               faceToNode(iface,1)=inode + 1;
00419               faceToNode(iface,2)=inode + (NX+1)*(NY+1)+1;
00420               faceToNode(iface,3)=inode + (NX+1)*(NY+1);
00421               if (j < NY) {
00422                  ielem=i+j*NX+k*NX*NY;
00423                  faceToEdge(iface,0)=elemToEdge(ielem,0);
00424                  faceToEdge(iface,1)=elemToEdge(ielem,9);
00425                  faceToEdge(iface,2)=elemToEdge(ielem,4);
00426                  faceToEdge(iface,3)=elemToEdge(ielem,8);
00427                  elemToFace(ielem,0)=iface;
00428                  if (j > 0) {
00429                     elemToFace(ielem-NX,2)=iface;
00430                  }
00431               }
00432               else if (j == NY) {
00433                  ielem=i+(NY-1)*NX+k*NX*NY;
00434                  faceToEdge(iface,0)=elemToEdge(ielem,2);
00435                  faceToEdge(iface,1)=elemToEdge(ielem,10);
00436                  faceToEdge(iface,2)=elemToEdge(ielem,6);
00437                  faceToEdge(iface,3)=elemToEdge(ielem,11);
00438                  elemToFace(ielem,2)=iface;
00439               }
00440               iface++;
00441            }
00442            if (j < NY && k < NZ) {
00443               faceToNode(iface,0)=inode;
00444               faceToNode(iface,1)=inode + NX+1;
00445               faceToNode(iface,2)=inode + (NX+1)*(NY+1) + NX+1;
00446               faceToNode(iface,3)=inode + (NX+1)*(NY+1);
00447               if (i < NX) {
00448                  ielem=i+j*NX+k*NX*NY;
00449                  faceToEdge(iface,0)=elemToEdge(ielem,3);
00450                  faceToEdge(iface,1)=elemToEdge(ielem,11);
00451                  faceToEdge(iface,2)=elemToEdge(ielem,7);
00452                  faceToEdge(iface,3)=elemToEdge(ielem,8);
00453                  elemToFace(ielem,3)=iface;
00454                  if (i > 0) {
00455                     elemToFace(ielem-1,1)=iface;
00456                  }
00457               }
00458               else if (i == NX) {
00459                  ielem=NX-1+j*NX+k*NX*NY;
00460                  faceToEdge(iface,0)=elemToEdge(ielem,1);
00461                  faceToEdge(iface,1)=elemToEdge(ielem,10);
00462                  faceToEdge(iface,2)=elemToEdge(ielem,5);
00463                  faceToEdge(iface,3)=elemToEdge(ielem,9);
00464                  elemToFace(ielem,1)=iface;
00465               }
00466               iface++;
00467            }
00468            if (i < NX && j < NY) {
00469               faceToNode(iface,0)=inode;
00470               faceToNode(iface,1)=inode + 1;
00471               faceToNode(iface,2)=inode + NX+2;
00472               faceToNode(iface,3)=inode + NX+1;
00473               if (k < NZ) {
00474                  ielem=i+j*NX+k*NX*NY;
00475                  faceToEdge(iface,0)=elemToEdge(ielem,0);
00476                  faceToEdge(iface,1)=elemToEdge(ielem,1);
00477                  faceToEdge(iface,2)=elemToEdge(ielem,2);
00478                  faceToEdge(iface,3)=elemToEdge(ielem,3);
00479                  elemToFace(ielem,4)=iface;
00480                  if (k > 0) {
00481                     elemToFace(ielem-NX*NY,5)=iface;
00482                  }
00483               }
00484               else if (k == NZ) {
00485                  ielem=i+j*NX+(NZ-1)*NX*NY;
00486                  faceToEdge(iface,0)=elemToEdge(ielem,4);
00487                  faceToEdge(iface,1)=elemToEdge(ielem,5);
00488                  faceToEdge(iface,2)=elemToEdge(ielem,6);
00489                  faceToEdge(iface,3)=elemToEdge(ielem,7);
00490                  elemToFace(ielem,5)=iface;
00491               }
00492               iface++;
00493            }
00494           inode++;
00495          }
00496       }
00497    }
00498 
00499    // Find boundary faces  
00500     FieldContainer<int> faceOnBoundary(numFaces);
00501     for (int i=0; i<numFaces; i++){
00502        if (nodeOnBoundary(faceToNode(i,0)) && nodeOnBoundary(faceToNode(i,1))
00503           && nodeOnBoundary(faceToNode(i,2)) && nodeOnBoundary(faceToNode(i,3))){
00504            faceOnBoundary(i)=1;
00505        }
00506     }
00507         
00508  
00509 #define DUMP_DATA
00510 #ifdef DUMP_DATA
00511    // Output connectivity
00512     ofstream fe2nout("elem2node.dat");
00513     ofstream fe2eout("elem2edge.dat");
00514     for (int k=0; k<NZ; k++) {
00515       for (int j=0; j<NY; j++) {
00516         for (int i=0; i<NX; i++) {
00517           int ielem = i + j * NX + k * NX * NY;
00518           for (int m=0; m<numNodesPerElem; m++){
00519               fe2nout << elemToNode(ielem,m) <<"  ";
00520            }
00521           fe2nout <<"\n";
00522           for (int l=0; l<numEdgesPerElem; l++) {
00523              fe2eout << elemToEdge(ielem,l) << "  ";
00524           }
00525           fe2eout << "\n";
00526         }
00527       }
00528     }
00529     fe2nout.close();
00530     fe2eout.close();
00531 #endif
00532 
00533 #ifdef DUMP_DATA_EXTRA
00534     ofstream fed2nout("edge2node.dat");
00535     for (int i=0; i<numEdges; i++) {
00536        fed2nout << edgeToNode(i,0) <<" ";
00537        fed2nout << edgeToNode(i,1) <<"\n";
00538     }
00539     fed2nout.close();
00540 
00541     ofstream fBnodeout("nodeOnBndy.dat");
00542     ofstream fBedgeout("edgeOnBndy.dat");
00543     for (int i=0; i<numNodes; i++) {
00544         fBnodeout << nodeOnBoundary(i) <<"\n";
00545     }
00546     for (int i=0; i<numEdges; i++) {
00547         fBedgeout << edgeOnBoundary(i) <<"\n";
00548     }
00549     fBnodeout.close();
00550     fBedgeout.close();
00551 #endif
00552 
00553    // Set material properties using undeformed grid assuming each element has only one value of mu
00554     FieldContainer<double> muVal(numElems);
00555     for (int k=0; k<NZ; k++) {
00556       for (int j=0; j<NY; j++) {
00557         for (int i=0; i<NX; i++) {
00558           int ielem = i + j * NX + k * NX * NY;
00559           double midElemX = nodeCoord(elemToNode(ielem,0),0) + hx/2.0;
00560           double midElemY = nodeCoord(elemToNode(ielem,0),1) + hy/2.0;
00561           double midElemZ = nodeCoord(elemToNode(ielem,0),2) + hz/2.0;
00562           if ( (midElemX > mu1LeftX) && (midElemY > mu1LeftY) && (midElemZ > mu1LeftZ) &&
00563                (midElemX <= mu1RightX) && (midElemY <= mu1RightY) && (midElemZ <= mu1RightZ) ){
00564              muVal(ielem) = mu1;
00565           }
00566            else {
00567              muVal(ielem) = mu2;
00568           }
00569         }
00570       }
00571     }
00572 
00573    // Perturb mesh coordinates (only interior nodes)
00574     if (randomMesh){
00575       for (int k=1; k<NZ; k++) {
00576         for (int j=1; j<NY; j++) {
00577           for (int i=1; i<NX; i++) {
00578             int inode = i + j * (NX + 1) + k * (NX + 1) * (NY + 1);
00579            // random numbers between -1.0 and 1.0
00580             double rx = 2.0 * (double)rand()/RAND_MAX - 1.0;
00581             double ry = 2.0 * (double)rand()/RAND_MAX - 1.0;
00582             double rz = 2.0 * (double)rand()/RAND_MAX - 1.0; 
00583            // limit variation to 1/4 edge length
00584             nodeCoord(inode,0) = nodeCoord(inode,0) + 0.125 * hx * rx;
00585             nodeCoord(inode,1) = nodeCoord(inode,1) + 0.125 * hy * ry;
00586             nodeCoord(inode,2) = nodeCoord(inode,2) + 0.125 * hz * rz;
00587           }
00588         }
00589       }
00590     }
00591 
00592 #ifdef DUMP_DATA
00593    // Print nodal coords
00594     ofstream fcoordout("coords.dat");
00595     for (int i=0; i<numNodes; i++) {
00596        fcoordout << nodeCoord(i,0) <<" ";
00597        fcoordout << nodeCoord(i,1) <<" ";
00598        fcoordout << nodeCoord(i,2) <<"\n";
00599     }
00600     fcoordout.close();
00601 #endif
00602 
00603 
00604 // **************************** INCIDENCE MATRIX **************************************
00605 
00606    // Node to edge incidence matrix
00607     *outStream << "Building incidence matrix ... \n\n";
00608 
00609     Epetra_SerialComm Comm;
00610     Epetra_Map globalMapC(numEdges, 0, Comm);
00611     Epetra_Map globalMapG(numNodes, 0, Comm);
00612     Epetra_FECrsMatrix DGrad(Copy, globalMapC, globalMapG, 2);
00613 
00614     double vals[2];
00615     vals[0]=-1.0; vals[1]=1.0;
00616     for (int j=0; j<numEdges; j++){
00617         int rowNum = j;
00618         int colNum[2];
00619         colNum[0] = edgeToNode(j,0);
00620         colNum[1] = edgeToNode(j,1);
00621         DGrad.InsertGlobalValues(1, &rowNum, 2, colNum, vals);
00622     }
00623 
00624 
00625 // ************************************ CUBATURE **************************************
00626 
00627    // Get numerical integration points and weights for element
00628     *outStream << "Getting cubature ... \n\n";
00629 
00630     DefaultCubatureFactory<double>  cubFactory;                                   
00631     int cubDegree = 2;
00632     Teuchos::RCP<Cubature<double> > hexCub = cubFactory.create(hex_8, cubDegree); 
00633 
00634     int cubDim       = hexCub->getDimension();
00635     int numCubPoints = hexCub->getNumPoints();
00636 
00637     FieldContainer<double> cubPoints(numCubPoints, cubDim);
00638     FieldContainer<double> cubWeights(numCubPoints);
00639 
00640     hexCub->getCubature(cubPoints, cubWeights);
00641 
00642 
00643    // Get numerical integration points and weights for hexahedron face
00644     //             (needed for rhs boundary term)
00645 
00646     // Define topology of the face parametrization domain as [-1,1]x[-1,1]
00647     CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00648 
00649     // Define cubature
00650     DefaultCubatureFactory<double>  cubFactoryFace;
00651     Teuchos::RCP<Cubature<double> > hexFaceCubature = cubFactoryFace.create(paramQuadFace, 3);
00652     int cubFaceDim    = hexFaceCubature -> getDimension();
00653     int numFacePoints = hexFaceCubature -> getNumPoints();
00654 
00655     // Define storage for cubature points and weights on [-1,1]x[-1,1]
00656     FieldContainer<double> paramGaussWeights(numFacePoints);
00657     FieldContainer<double> paramGaussPoints(numFacePoints,cubFaceDim);
00658 
00659     // Define storage for cubature points on workset faces
00660     hexFaceCubature -> getCubature(paramGaussPoints, paramGaussWeights);
00661 
00662 
00663 // ************************************** BASIS ***************************************
00664 
00665    // Define basis 
00666     *outStream << "Getting basis ... \n\n";
00667     Basis_HCURL_HEX_I1_FEM<double, FieldContainer<double> > hexHCurlBasis;
00668     Basis_HGRAD_HEX_C1_FEM<double, FieldContainer<double> > hexHGradBasis;
00669 
00670     int numFieldsC = hexHCurlBasis.getCardinality();
00671     int numFieldsG = hexHGradBasis.getCardinality();
00672 
00673   // Evaluate basis at cubature points
00674      FieldContainer<double> hexGVals(numFieldsG, numCubPoints); 
00675      FieldContainer<double> hexCVals(numFieldsC, numCubPoints, spaceDim); 
00676      FieldContainer<double> hexCurls(numFieldsC, numCubPoints, spaceDim); 
00677      FieldContainer<double> worksetCVals(numFieldsC, numFacePoints, spaceDim);
00678 
00679 
00680      hexHCurlBasis.getValues(hexCVals, cubPoints, OPERATOR_VALUE);
00681      hexHCurlBasis.getValues(hexCurls, cubPoints, OPERATOR_CURL);
00682      hexHGradBasis.getValues(hexGVals, cubPoints, OPERATOR_VALUE);
00683 
00684 
00685 // ******** LOOP OVER ELEMENTS TO CREATE LOCAL MASS and STIFFNESS MATRICES *************
00686 
00687 
00688     *outStream << "Building mass and stiffness matrices ... \n\n";
00689 
00690  // Settings and data structures for mass and stiffness matrices
00691     typedef CellTools<double>  CellTools;
00692     typedef FunctionSpaceTools fst;
00693     typedef ArrayTools art;
00694     int numCells = 1; 
00695 
00696    // Containers for nodes and edge signs 
00697     FieldContainer<double> hexNodes(numCells, numNodesPerElem, spaceDim);
00698     FieldContainer<double> hexEdgeSigns(numCells, numFieldsC);
00699    // Containers for Jacobian
00700     FieldContainer<double> hexJacobian(numCells, numCubPoints, spaceDim, spaceDim);
00701     FieldContainer<double> hexJacobInv(numCells, numCubPoints, spaceDim, spaceDim);
00702     FieldContainer<double> hexJacobDet(numCells, numCubPoints);
00703    // Containers for element HGRAD mass matrix
00704     FieldContainer<double> massMatrixG(numCells, numFieldsG, numFieldsG);
00705     FieldContainer<double> weightedMeasure(numCells, numCubPoints);
00706     FieldContainer<double> weightedMeasureMuInv(numCells, numCubPoints);
00707     FieldContainer<double> hexGValsTransformed(numCells, numFieldsG, numCubPoints);
00708     FieldContainer<double> hexGValsTransformedWeighted(numCells, numFieldsG, numCubPoints);
00709    // Containers for element HCURL mass matrix
00710     FieldContainer<double> massMatrixC(numCells, numFieldsC, numFieldsC);
00711     FieldContainer<double> hexCValsTransformed(numCells, numFieldsC, numCubPoints, spaceDim);
00712     FieldContainer<double> hexCValsTransformedWeighted(numCells, numFieldsC, numCubPoints, spaceDim);
00713    // Containers for element HCURL stiffness matrix
00714     FieldContainer<double> stiffMatrixC(numCells, numFieldsC, numFieldsC);
00715     FieldContainer<double> weightedMeasureMu(numCells, numCubPoints);    
00716     FieldContainer<double> hexCurlsTransformed(numCells, numFieldsC, numCubPoints, spaceDim);
00717     FieldContainer<double> hexCurlsTransformedWeighted(numCells, numFieldsC, numCubPoints, spaceDim);
00718    // Containers for right hand side vectors
00719     FieldContainer<double> rhsDatag(numCells, numCubPoints, cubDim);
00720     FieldContainer<double> rhsDatah(numCells, numCubPoints, cubDim);
00721     FieldContainer<double> gC(numCells, numFieldsC);
00722     FieldContainer<double> hC(numCells, numFieldsC);
00723     FieldContainer<double> hCBoundary(numCells, numFieldsC);
00724     FieldContainer<double> refGaussPoints(numFacePoints,spaceDim);
00725     FieldContainer<double> worksetGaussPoints(numCells,numFacePoints,spaceDim);
00726     FieldContainer<double> worksetJacobians(numCells, numFacePoints, spaceDim, spaceDim);
00727     FieldContainer<double> worksetJacobInv(numCells, numFacePoints, spaceDim, spaceDim);
00728     FieldContainer<double> worksetFaceN(numCells, numFacePoints, spaceDim);
00729     FieldContainer<double> worksetFaceNweighted(numCells, numFacePoints, spaceDim);
00730     FieldContainer<double> worksetVFieldVals(numCells, numFacePoints, spaceDim);
00731     FieldContainer<double> worksetCValsTransformed(numCells, numFieldsC, numFacePoints, spaceDim);
00732     FieldContainer<double> divuFace(numCells, numFacePoints);
00733     FieldContainer<double> worksetFieldDotNormal(numCells, numFieldsC, numFacePoints);
00734    // Container for cubature points in physical space
00735     FieldContainer<double> physCubPoints(numCells,numCubPoints, cubDim);
00736 
00737     
00738    // Global arrays in Epetra format
00739     Epetra_FECrsMatrix MassG(Copy, globalMapG, numFieldsG);
00740     Epetra_FECrsMatrix MassC(Copy, globalMapC, numFieldsC);
00741     Epetra_FECrsMatrix StiffC(Copy, globalMapC, numFieldsC);
00742     Epetra_FEVector rhsC(globalMapC);
00743 
00744 #ifdef DUMP_DATA
00745     ofstream fSignsout("edgeSigns.dat");
00746 #endif
00747 
00748  // *** Element loop ***
00749     for (int k=0; k<numElems; k++) {
00750 
00751      // Physical cell coordinates
00752       for (int i=0; i<numNodesPerElem; i++) {
00753          hexNodes(0,i,0) = nodeCoord(elemToNode(k,i),0);
00754          hexNodes(0,i,1) = nodeCoord(elemToNode(k,i),1);
00755          hexNodes(0,i,2) = nodeCoord(elemToNode(k,i),2);
00756       }
00757 
00758      // Edge signs
00759       for (int j=0; j<numEdgesPerElem; j++) {
00760           if (elemToNode(k,refEdgeToNode(j,0))==edgeToNode(elemToEdge(k,j),0) &&
00761               elemToNode(k,refEdgeToNode(j,1))==edgeToNode(elemToEdge(k,j),1))
00762               hexEdgeSigns(0,j) = 1.0;
00763           else 
00764               hexEdgeSigns(0,j) = -1.0;
00765 #ifdef DUMP_DATA
00766           fSignsout << hexEdgeSigns(0,j) << "  ";
00767 #endif
00768       }
00769 #ifdef DUMP_DATA
00770        fSignsout << "\n";
00771 #endif
00772 
00773     // Compute cell Jacobians, their inverses and their determinants
00774        CellTools::setJacobian(hexJacobian, cubPoints, hexNodes, hex_8);
00775        CellTools::setJacobianInv(hexJacobInv, hexJacobian );
00776        CellTools::setJacobianDet(hexJacobDet, hexJacobian );
00777 
00778 // ************************** Compute element HGrad mass matrices *******************************
00779   
00780      // transform to physical coordinates 
00781       fst::HGRADtransformVALUE<double>(hexGValsTransformed, hexGVals);
00782       
00783      // compute weighted measure
00784       fst::computeCellMeasure<double>(weightedMeasure, hexJacobDet, cubWeights);
00785 
00786       // combine mu value with weighted measure
00787       for (int nC = 0; nC < numCells; nC++){
00788         for (int nPt = 0; nPt < numCubPoints; nPt++){
00789           weightedMeasureMuInv(nC,nPt) = weightedMeasure(nC,nPt) / muVal(k);
00790         }
00791       }
00792       
00793      // multiply values with weighted measure
00794       fst::multiplyMeasure<double>(hexGValsTransformedWeighted,
00795                                    weightedMeasureMuInv, hexGValsTransformed);
00796 
00797      // integrate to compute element mass matrix
00798       fst::integrate<double>(massMatrixG,
00799                              hexGValsTransformed, hexGValsTransformedWeighted, COMP_CPP);
00800 
00801       // assemble into global matrix
00802       int err = 0;
00803       for (int row = 0; row < numFieldsG; row++){
00804         for (int col = 0; col < numFieldsG; col++){
00805             int rowIndex = elemToNode(k,row);
00806             int colIndex = elemToNode(k,col);
00807             double val = massMatrixG(0,row,col);
00808             MassG.InsertGlobalValues(1, &rowIndex, 1, &colIndex, &val);
00809          }
00810       }
00811 
00812 // ************************** Compute element HCurl mass matrices *******************************
00813 
00814      // transform to physical coordinates 
00815       fst::HCURLtransformVALUE<double>(hexCValsTransformed, hexJacobInv, 
00816                                    hexCVals);
00817 
00818      // multiply by weighted measure
00819       fst::multiplyMeasure<double>(hexCValsTransformedWeighted,
00820                                    weightedMeasure, hexCValsTransformed);
00821 
00822      // integrate to compute element mass matrix
00823       fst::integrate<double>(massMatrixC,
00824                              hexCValsTransformed, hexCValsTransformedWeighted,
00825                              COMP_CPP);
00826 
00827      // apply edge signs
00828       fst::applyLeftFieldSigns<double>(massMatrixC, hexEdgeSigns);
00829       fst::applyRightFieldSigns<double>(massMatrixC, hexEdgeSigns);
00830 
00831 
00832      // assemble into global matrix
00833       err = 0;
00834       for (int row = 0; row < numFieldsC; row++){
00835         for (int col = 0; col < numFieldsC; col++){
00836             int rowIndex = elemToEdge(k,row);
00837             int colIndex = elemToEdge(k,col);
00838             double val = massMatrixC(0,row,col);
00839             MassC.InsertGlobalValues(1, &rowIndex, 1, &colIndex, &val);
00840          }
00841       }
00842 
00843 // ************************ Compute element HCurl stiffness matrices *****************************
00844 
00845       // transform to physical coordinates 
00846       fst::HCURLtransformCURL<double>(hexCurlsTransformed, hexJacobian, hexJacobDet, 
00847                                        hexCurls);
00848 
00849       // combine mu value with weighted measure
00850       for (int nC = 0; nC < numCells; nC++){
00851         for (int nPt = 0; nPt < numCubPoints; nPt++){
00852           weightedMeasureMu(nC,nPt) = weightedMeasure(nC,nPt) / muVal(k);
00853          }
00854       }
00855 
00856      // multiply by weighted measure
00857       fst::multiplyMeasure<double>(hexCurlsTransformedWeighted,
00858                                    weightedMeasureMu, hexCurlsTransformed);
00859 
00860      // integrate to compute element stiffness matrix
00861       fst::integrate<double>(stiffMatrixC,
00862                              hexCurlsTransformed, hexCurlsTransformedWeighted,
00863                              COMP_CPP);
00864 
00865      // apply edge signs
00866      fst::applyLeftFieldSigns<double>(stiffMatrixC, hexEdgeSigns);
00867      fst::applyRightFieldSigns<double>(stiffMatrixC, hexEdgeSigns);
00868 
00869      // assemble into global matrix
00870       err = 0;
00871       for (int row = 0; row < numFieldsC; row++){
00872         for (int col = 0; col < numFieldsC; col++){
00873             int rowIndex = elemToEdge(k,row);
00874             int colIndex = elemToEdge(k,col);
00875             double val = stiffMatrixC(0,row,col);
00876             StiffC.InsertGlobalValues(1, &rowIndex, 1, &colIndex, &val);
00877          }
00878       }
00879 
00880 // ******************************* Build right hand side ************************************
00881 
00882       // transform integration points to physical points
00883        FieldContainer<double> physCubPoints(numCells,numCubPoints, cubDim);
00884        CellTools::mapToPhysicalFrame(physCubPoints, cubPoints, hexNodes, hex_8);
00885 
00886       // evaluate right hand side functions at physical points
00887        FieldContainer<double> rhsDatag(numCells, numCubPoints, cubDim);
00888        FieldContainer<double> rhsDatah(numCells, numCubPoints, cubDim);
00889        for (int nPt = 0; nPt < numCubPoints; nPt++){
00890 
00891           double x = physCubPoints(0,nPt,0);
00892           double y = physCubPoints(0,nPt,1);
00893           double z = physCubPoints(0,nPt,2);
00894           double du1, du2, du3;
00895 
00896           evalCurlu(du1, du2, du3, x, y, z);
00897           rhsDatag(0,nPt,0) = du1;
00898           rhsDatag(0,nPt,1) = du2;
00899           rhsDatag(0,nPt,2) = du3;
00900          
00901           evalGradDivu(du1, du2, du3,  x, y, z);
00902           rhsDatah(0,nPt,0) = du1;
00903           rhsDatah(0,nPt,1) = du2;
00904           rhsDatah(0,nPt,2) = du3;
00905        }
00906 
00907      // integrate (g,curl w) term
00908       fst::integrate<double>(gC, rhsDatag, hexCurlsTransformedWeighted,
00909                              COMP_CPP);
00910 
00911      // integrate -(grad h, w)  
00912       fst::integrate<double>(hC, rhsDatah, hexCValsTransformedWeighted,
00913                              COMP_CPP);
00914 
00915      // apply signs
00916       fst::applyFieldSigns<double>(gC, hexEdgeSigns);
00917       fst::applyFieldSigns<double>(hC, hexEdgeSigns);
00918 
00919 
00920      // calculate boundary term, (h*w, n)_{\Gamma} 
00921       for (int i = 0; i < numFacesPerElem; i++){
00922         if (faceOnBoundary(elemToFace(k,i))){
00923 
00924          // Map Gauss points on quad to reference face: paramGaussPoints -> refGaussPoints
00925             CellTools::mapToReferenceSubcell(refGaussPoints,
00926                                    paramGaussPoints,
00927                                    2, i, hex_8);
00928 
00929          // Get basis values at points on reference cell
00930            hexHCurlBasis.getValues(worksetCVals, refGaussPoints, OPERATOR_VALUE);
00931 
00932          // Compute Jacobians at Gauss pts. on reference face for all parent cells
00933            CellTools::setJacobian(worksetJacobians,
00934                          refGaussPoints,
00935                          hexNodes, hex_8);
00936            CellTools::setJacobianInv(worksetJacobInv, worksetJacobians );
00937 
00938          // transform to physical coordinates
00939             fst::HCURLtransformVALUE<double>(worksetCValsTransformed, worksetJacobInv,
00940                                    worksetCVals);
00941 
00942          // Map Gauss points on quad from ref. face to face workset: refGaussPoints -> worksetGaussPoints
00943             CellTools::mapToPhysicalFrame(worksetGaussPoints,
00944                                 refGaussPoints,
00945                                 hexNodes, hex_8);
00946 
00947          // Compute face normals
00948             CellTools::getPhysicalFaceNormals(worksetFaceN,
00949                                               worksetJacobians,
00950                                               i, hex_8);
00951 
00952          // multiply with weights
00953             for(int nPt = 0; nPt < numFacePoints; nPt++){
00954                 for (int dim = 0; dim < spaceDim; dim++){
00955                    worksetFaceNweighted(0,nPt,dim) = worksetFaceN(0,nPt,dim) * paramGaussWeights(nPt);
00956                 } //dim
00957              } //nPt
00958 
00959             fst::dotMultiplyDataField<double>(worksetFieldDotNormal, worksetFaceNweighted, 
00960                                                worksetCValsTransformed);
00961 
00962          // Evaluate div u at face points
00963            for(int nPt = 0; nPt < numFacePoints; nPt++){
00964 
00965              double x = worksetGaussPoints(0, nPt, 0);
00966              double y = worksetGaussPoints(0, nPt, 1);
00967              double z = worksetGaussPoints(0, nPt, 2);
00968 
00969              divuFace(0,nPt)=evalDivu(x, y, z);
00970            }
00971 
00972           // Integrate
00973           fst::integrate<double>(hCBoundary, divuFace, worksetFieldDotNormal,
00974                              COMP_CPP);
00975 
00976           // apply signs
00977            fst::applyFieldSigns<double>(hCBoundary, hexEdgeSigns);
00978 
00979          // add into hC term
00980             for (int nF = 0; nF < numFieldsC; nF++){
00981                 hC(0,nF) = hC(0,nF) - hCBoundary(0,nF);
00982             }
00983 
00984         } // if faceOnBoundary
00985       } // numFaces
00986 
00987 
00988     // assemble into global vector
00989      for (int row = 0; row < numFieldsC; row++){
00990            int rowIndex = elemToEdge(k,row);
00991            double val = gC(0,row)-hC(0,row);
00992            rhsC.SumIntoGlobalValues(1, &rowIndex, &val);
00993      }
00994  
00995      
00996  } // *** end element loop ***
00997 
00998 #ifdef DUMP_DATA
00999    fSignsout.close();
01000 #endif
01001 
01002   // Assemble over multiple processors, if necessary
01003    MassG.GlobalAssemble();  MassG.FillComplete();
01004    MassC.GlobalAssemble();  MassC.FillComplete();
01005    StiffC.GlobalAssemble(); StiffC.FillComplete();
01006    rhsC.GlobalAssemble();
01007    DGrad.GlobalAssemble(); DGrad.FillComplete(MassG.RowMap(),MassC.RowMap());
01008 
01009 
01010   // Build the inverse diagonal for MassG
01011    Epetra_CrsMatrix MassGinv(Copy,MassG.RowMap(),MassG.RowMap(),1);
01012    Epetra_Vector DiagG(MassG.RowMap());
01013 
01014    DiagG.PutScalar(1.0);
01015    MassG.Multiply(false,DiagG,DiagG);
01016    for(int i=0; i<DiagG.MyLength(); i++) {
01017      DiagG[i]=1.0/DiagG[i];
01018    }
01019    for(int i=0; i<DiagG.MyLength(); i++) {
01020      int CID=MassG.GCID(i);
01021      MassGinv.InsertGlobalValues(MassG.GRID(i),1,&(DiagG[i]),&CID);
01022    }
01023    MassGinv.FillComplete();
01024 
01025   // Set value to zero on diagonal that cooresponds to boundary node
01026    for(int i=0;i<numNodes;i++) {
01027      if (nodeOnBoundary(i)){
01028       double val=0.0;
01029       MassGinv.ReplaceGlobalValues(i,1,&val,&i);
01030      }
01031    }
01032 
01033     int numEntries;
01034     double *values;
01035     int *cols;
01036   
01037   // Adjust matrices and rhs due to boundary conditions
01038    for (int row = 0; row<numEdges; row++){
01039       MassC.ExtractMyRowView(row,numEntries,values,cols);
01040         for (int i=0; i<numEntries; i++){
01041            if (edgeOnBoundary(cols[i])) {
01042              values[i]=0;
01043           }
01044        }
01045       StiffC.ExtractMyRowView(row,numEntries,values,cols);
01046         for (int i=0; i<numEntries; i++){
01047            if (edgeOnBoundary(cols[i])) {
01048              values[i]=0;
01049           }
01050        }
01051     }
01052    for (int row = 0; row<numEdges; row++){
01053        if (edgeOnBoundary(row)) {
01054           int rowindex = row;
01055           StiffC.ExtractMyRowView(row,numEntries,values,cols);
01056           for (int i=0; i<numEntries; i++){
01057              values[i]=0;
01058           }
01059           MassC.ExtractMyRowView(row,numEntries,values,cols);
01060           for (int i=0; i<numEntries; i++){
01061              values[i]=0;
01062           }
01063          rhsC[0][row]=0;
01064          double val = 1.0;
01065          StiffC.ReplaceGlobalValues(1, &rowindex, 1, &rowindex, &val);
01066        }
01067     }
01068 
01069 
01070 #ifdef DUMP_DATA
01071   // Dump matrices to disk
01072    EpetraExt::RowMatrixToMatlabFile("mag_m0inv_matrix.dat",MassGinv);
01073    EpetraExt::RowMatrixToMatlabFile("mag_m1_matrix.dat",MassC);
01074    EpetraExt::RowMatrixToMatlabFile("mag_k1_matrix.dat",StiffC);
01075    EpetraExt::RowMatrixToMatlabFile("mag_t0_matrix.dat",DGrad);
01076    EpetraExt::MultiVectorToMatrixMarketFile("mag_rhs1_vector.dat",rhsC,0,0,false);
01077    fSignsout.close();
01078 #endif
01079 
01080 
01081    std::cout << "End Result: TEST PASSED\n";
01082 
01083  // reset format state of std::cout
01084  std::cout.copyfmt(oldFormatState);
01085  
01086  return 0;
01087 }
01088 
01089 // Calculates value of exact solution u
01090  int evalu(double & uExact0, double & uExact1, double & uExact2, double & x, double & y, double & z)
01091  {
01092    // function 1
01093     uExact0 = exp(x+y+z)*(y+1.0)*(y-1.0)*(z+1.0)*(z-1.0);
01094     uExact1 = exp(x+y+z)*(x+1.0)*(x-1.0)*(z+1.0)*(z-1.0);
01095     uExact2 = exp(x+y+z)*(x+1.0)*(x-1.0)*(y+1.0)*(y-1.0);
01096 
01097  /*
01098    // function 2
01099     uExact0 = cos(M_PI*x)*exp(y*z)*(y+1.0)*(y-1.0)*(z+1.0)*(z-1.0);
01100     uExact1 = cos(M_PI*y)*exp(x*z)*(x+1.0)*(x-1.0)*(z+1.0)*(z-1.0);
01101     uExact2 = cos(M_PI*z)*exp(x*y)*(x+1.0)*(x-1.0)*(y+1.0)*(y-1.0);
01102 
01103    // function 3
01104     uExact0 = cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
01105     uExact1 = sin(M_PI*x)*cos(M_PI*y)*sin(M_PI*z);
01106     uExact2 = sin(M_PI*x)*sin(M_PI*y)*cos(M_PI*z);
01107  
01108    // function 4
01109     uExact0 = (y*y - 1.0)*(z*z-1.0);
01110     uExact1 = (x*x - 1.0)*(z*z-1.0);
01111     uExact2 = (x*x - 1.0)*(y*y-1.0);
01112  */
01113   
01114    return 0;
01115  }
01116 
01117 // Calculates divergence of exact solution u
01118  double evalDivu(double & x, double & y, double & z)
01119  {
01120    // function 1
01121     double divu = exp(x+y+z)*(y+1.0)*(y-1.0)*(z+1.0)*(z-1.0)
01122                  + exp(x+y+z)*(x+1.0)*(x-1.0)*(z+1.0)*(z-1.0)
01123                  + exp(x+y+z)*(x+1.0)*(x-1.0)*(y+1.0)*(y-1.0);
01124 
01125  /*
01126    // function 2
01127     double divu = -M_PI*sin(M_PI*x)*exp(y*z)*(y+1.0)*(y-1.0)*(z+1.0)*(z-1.0)
01128                  -M_PI*sin(M_PI*y)*exp(x*z)*(x+1.0)*(x-1.0)*(z+1.0)*(z-1.0)
01129                  -M_PI*sin(M_PI*z)*exp(x*y)*(x+1.0)*(x-1.0)*(y+1.0)*(y-1.0);
01130 
01131    // function 3
01132     double divu = -3.0*M_PI*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
01133 
01134    // function 4
01135     double divu = 0.0;
01136   */
01137 
01138    return divu;
01139  }
01140 
01141 
01142 // Calculates curl of exact solution u
01143  int evalCurlu(double & curlu0, double & curlu1, double & curlu2, double & x, double & y, double & z)
01144  {
01145   
01146    // function 1
01147     double duxdy = exp(x+y+z)*(z*z-1.0)*(y*y+2.0*y-1.0);
01148     double duxdz = exp(x+y+z)*(y*y-1.0)*(z*z+2.0*z-1.0);
01149     double duydx = exp(x+y+z)*(z*z-1.0)*(x*x+2.0*x-1.0);
01150     double duydz = exp(x+y+z)*(x*x-1.0)*(z*z+2.0*z-1.0);
01151     double duzdx = exp(x+y+z)*(y*y-1.0)*(x*x+2.0*x-1.0);
01152     double duzdy = exp(x+y+z)*(x*x-1.0)*(y*y+2.0*y-1.0);
01153 
01154  /*
01155    // function 2
01156     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);
01157     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);
01158     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);
01159     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);
01160     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);
01161     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);
01162 
01163    // function 3
01164     double duxdy = M_PI*cos(M_PI*x)*cos(M_PI*y)*sin(M_PI*z);
01165     double duxdz = M_PI*cos(M_PI*x)*sin(M_PI*y)*cos(M_PI*z);
01166     double duydx = M_PI*cos(M_PI*x)*cos(M_PI*y)*sin(M_PI*z);
01167     double duydz = M_PI*sin(M_PI*x)*cos(M_PI*y)*cos(M_PI*z);
01168     double duzdx = M_PI*cos(M_PI*x)*sin(M_PI*y)*cos(M_PI*z);
01169     double duzdy = M_PI*sin(M_PI*x)*cos(M_PI*y)*cos(M_PI*z);
01170  
01171    // function 4
01172     double duxdy = 2.0*y*(z*z-1);
01173     double duxdz = 2.0*z*(y*y-1);
01174     double duydx = 2.0*x*(z*z-1);
01175     double duydz = 2.0*z*(x*x-1);
01176     double duzdx = 2.0*x*(y*y-1);
01177     double duzdy = 2.0*y*(x*x-1);
01178 */
01179 
01180    curlu0 = duzdy - duydz;
01181    curlu1 = duxdz - duzdx;
01182    curlu2 = duydx - duxdy;
01183 
01184    return 0;
01185  }
01186 
01187 // Calculates gradient of the divergence of exact solution u
01188  int evalGradDivu(double & gradDivu0, double & gradDivu1, double & gradDivu2, double & x, double & y, double & z)
01189 {
01190    
01191    // function 1
01192     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));
01193     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));
01194     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));
01195  
01196   /*
01197    // function 2
01198     gradDivu0 = -M_PI*M_PI*cos(M_PI*x)*exp(y*z)*(y+1.0)*(y-1.0)*(z+1.0)*(z-1.0)
01199                   -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)
01200                   -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);
01201     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)
01202                   -M_PI*M_PI*cos(M_PI*y)*exp(x*z)*(x+1.0)*(x-1.0)*(z+1.0)*(z-1.0)
01203                   -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);
01204     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)
01205                   -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)
01206                   -M_PI*M_PI*cos(M_PI*z)*exp(x*y)*(x+1.0)*(x-1.0)*(y+1.0)*(y-1.0);
01207 
01208    // function 3
01209     gradDivu0 = -3.0*M_PI*M_PI*cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
01210     gradDivu1 = -3.0*M_PI*M_PI*sin(M_PI*x)*cos(M_PI*y)*sin(M_PI*z);
01211     gradDivu2 = -3.0*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y)*cos(M_PI*z);
01212 
01213    // function 4
01214     gradDivu0 = 0;
01215     gradDivu1 = 0;
01216     gradDivu2 = 0;
01217  */
01218    
01219    return 0;
01220 }