|
Intrepid
|
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 }
1.7.6.1