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