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