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