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