Intrepid
example_02.cpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
93// Intrepid includes
99#include "Intrepid_HGRAD_HEX_C1_FEM.hpp"
103#include "Intrepid_Utils.hpp"
104
105// Epetra includes
106#include "Epetra_Time.h"
107#include "Epetra_Map.h"
108#include "Epetra_SerialComm.h"
109#include "Epetra_FECrsMatrix.h"
110#include "Epetra_FEVector.h"
111#include "Epetra_Vector.h"
112
113// Teuchos includes
114#include "Teuchos_oblackholestream.hpp"
115#include "Teuchos_RCP.hpp"
116#include "Teuchos_BLAS.hpp"
117
118// Shards includes
119#include "Shards_CellTopology.hpp"
120
121// EpetraExt includes
122#include "EpetraExt_RowMatrixOut.h"
123#include "EpetraExt_MultiVectorOut.h"
124
125using namespace std;
126using namespace Intrepid;
127
128// Functions to evaluate exact solution and its derivatives
129int evalu(double & uExact0, double & uExact1, double & uExact2, double & x, double & y, double & z);
130double evalDivu(double & x, double & y, double & z);
131int evalCurlu(double & curlu0, double & curlu1, double & curlu2, double & x, double & y, double & z);
132int evalCurlCurlu(double & curlCurlu0, double & curlCurlu1, double & curlCurlu2, double & x, double & y, double & z);
133
134int main(int argc, char *argv[]) {
135
136 //Check number of arguments
137 if (argc < 13) {
138 std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
139 std::cout <<"Usage:\n\n";
140 std::cout <<" ./Intrepid_example_Drivers_Example_02.exe NX NY NZ randomMesh mu1 mu2 mu1LX mu1RX mu1LY mu1RY mu1LZ mu1RZ verbose\n\n";
141 std::cout <<" where \n";
142 std::cout <<" int NX - num intervals in x direction (assumed box domain, -1,1) \n";
143 std::cout <<" int NY - num intervals in y direction (assumed box domain, -1,1) \n";
144 std::cout <<" int NZ - num intervals in z direction (assumed box domain, -1,1) \n";
145 std::cout <<" int randomMesh - 1 if mesh randomizer is to be used 0 if not \n";
146 std::cout <<" double mu1 - material property value for region 1 \n";
147 std::cout <<" double mu2 - material property value for region 2 \n";
148 std::cout <<" double mu1LX - left X boundary for region 1 \n";
149 std::cout <<" double mu1RX - right X boundary for region 1 \n";
150 std::cout <<" double mu1LY - left Y boundary for region 1 \n";
151 std::cout <<" double mu1RY - right Y boundary for region 1 \n";
152 std::cout <<" double mu1LZ - bottom Z boundary for region 1 \n";
153 std::cout <<" double mu1RZ - top Z boundary for region 1 \n";
154 std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
155 exit(1);
156 }
157
158 // This little trick lets us print to std::cout only if
159 // a (dummy) command-line argument is provided.
160 int iprint = argc - 1;
161 Teuchos::RCP<std::ostream> outStream;
162 Teuchos::oblackholestream bhs; // outputs nothing
163 if (iprint > 12)
164 outStream = Teuchos::rcp(&std::cout, false);
165 else
166 outStream = Teuchos::rcp(&bhs, false);
167
168 // Save the format state of the original std::cout.
169 Teuchos::oblackholestream oldFormatState;
170 oldFormatState.copyfmt(std::cout);
171
172 *outStream \
173 << "===============================================================================\n" \
174 << "| |\n" \
175 << "| Example: Generate Mass and Stiffness Matrices and Right-Hand Side Vector |\n"
176 << "| for Div-Curl System on Hexahedral Mesh with Div-Conforming Elements |\n" \
177 << "| |\n" \
178 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
179 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
180 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
181 << "| |\n" \
182 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
183 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
184 << "| |\n" \
185 << "===============================================================================\n";
186
187
188// ************************************ GET INPUTS **************************************
189
190 /* In the implementation for discontinuous material properties only the boundaries for
191 region 1, associated with mu1, are input. The remainder of the grid is assumed to use mu2.
192 Note that the material properties are assigned using the undeformed grid. */
193
194 int NX = atoi(argv[1]); // num intervals in x direction (assumed box domain, -1,1)
195 int NY = atoi(argv[2]); // num intervals in y direction (assumed box domain, -1,1)
196 int NZ = atoi(argv[3]); // num intervals in z direction (assumed box domain, -1,1)
197 int randomMesh = atoi(argv[4]); // 1 if mesh randomizer is to be used 0 if not
198 double mu1 = atof(argv[5]); // material property value for region 1
199 double mu2 = atof(argv[6]); // material property value for region 2
200 double mu1LeftX = atof(argv[7]); // left X boundary for region 1
201 double mu1RightX = atof(argv[8]); // right X boundary for region 1
202 double mu1LeftY = atof(argv[9]); // left Y boundary for region 1
203 double mu1RightY = atof(argv[10]); // right Y boundary for region 1
204 double mu1LeftZ = atof(argv[11]); // left Z boundary for region 1
205 double mu1RightZ = atof(argv[12]); // right Z boundary for region 1
206
207// *********************************** CELL TOPOLOGY **********************************
208
209 // Get cell topology for base hexahedron
210 typedef shards::CellTopology CellTopology;
211 CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() );
212
213 // Get dimensions
214 int numNodesPerElem = hex_8.getNodeCount();
215 int numEdgesPerElem = hex_8.getEdgeCount();
216 int numFacesPerElem = hex_8.getSideCount();
217 int numNodesPerEdge = 2;
218 int numNodesPerFace = 4;
219 int numEdgesPerFace = 4;
220 int spaceDim = hex_8.getDimension();
221
222 // Build reference element edge to node map
223 FieldContainer<int> refEdgeToNode(numEdgesPerElem,numNodesPerEdge);
224 for (int i=0; i<numEdgesPerElem; i++){
225 refEdgeToNode(i,0)=hex_8.getNodeMap(1, i, 0);
226 refEdgeToNode(i,1)=hex_8.getNodeMap(1, i, 1);
227 }
228
229 // Build reference element face to node map
230 FieldContainer<int> refFaceToNode(numFacesPerElem,numNodesPerFace);
231 for (int i=0; i<numFacesPerElem; i++){
232 refFaceToNode(i,0)=hex_8.getNodeMap(2, i, 0);
233 refFaceToNode(i,1)=hex_8.getNodeMap(2, i, 1);
234 refFaceToNode(i,2)=hex_8.getNodeMap(2, i, 2);
235 refFaceToNode(i,3)=hex_8.getNodeMap(2, i, 3);
236 }
237
238// *********************************** GENERATE MESH ************************************
239
240 *outStream << "Generating mesh ... \n\n";
241
242 *outStream << " NX" << " NY" << " NZ\n";
243 *outStream << std::setw(5) << NX <<
244 std::setw(5) << NY <<
245 std::setw(5) << NZ << "\n\n";
246
247 // Print mesh information
248 int numElems = NX*NY*NZ;
249 int numNodes = (NX+1)*(NY+1)*(NZ+1);
250 int numEdges = (NX)*(NY + 1)*(NZ + 1) + (NX + 1)*(NY)*(NZ + 1) + (NX + 1)*(NY + 1)*(NZ);
251 int numFaces = (NX)*(NY)*(NZ + 1) + (NX)*(NY + 1)*(NZ) + (NX + 1)*(NY)*(NZ);
252 *outStream << " Number of Elements: " << numElems << " \n";
253 *outStream << " Number of Nodes: " << numNodes << " \n";
254 *outStream << " Number of Edges: " << numEdges << " \n";
255 *outStream << " Number of Faces: " << numFaces << " \n\n";
256
257 // Cube
258 double leftX = -1.0, rightX = 1.0;
259 double leftY = -1.0, rightY = 1.0;
260 double leftZ = -1.0, rightZ = 1.0;
261
262 // Mesh spacing
263 double hx = (rightX-leftX)/((double)NX);
264 double hy = (rightY-leftY)/((double)NY);
265 double hz = (rightZ-leftZ)/((double)NZ);
266
267 // Get nodal coordinates
268 FieldContainer<double> nodeCoord(numNodes, spaceDim);
269 FieldContainer<int> nodeOnBoundary(numNodes);
270 int inode = 0;
271 for (int k=0; k<NZ+1; k++) {
272 for (int j=0; j<NY+1; j++) {
273 for (int i=0; i<NX+1; i++) {
274 nodeCoord(inode,0) = leftX + (double)i*hx;
275 nodeCoord(inode,1) = leftY + (double)j*hy;
276 nodeCoord(inode,2) = leftZ + (double)k*hz;
277 if (k==0 || j==0 || i==0 || k==NZ || j==NY || i==NX){
278 nodeOnBoundary(inode)=1;
279 }
280 inode++;
281 }
282 }
283 }
284
285 // Element to Node map
286 FieldContainer<int> elemToNode(numElems, numNodesPerElem);
287 int ielem = 0;
288 for (int k=0; k<NZ; k++) {
289 for (int j=0; j<NY; j++) {
290 for (int i=0; i<NX; i++) {
291 elemToNode(ielem,0) = (NY + 1)*(NX + 1)*k + (NX + 1)*j + i;
292 elemToNode(ielem,1) = (NY + 1)*(NX + 1)*k + (NX + 1)*j + i + 1;
293 elemToNode(ielem,2) = (NY + 1)*(NX + 1)*k + (NX + 1)*(j + 1) + i + 1;
294 elemToNode(ielem,3) = (NY + 1)*(NX + 1)*k + (NX + 1)*(j + 1) + i;
295 elemToNode(ielem,4) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*j + i;
296 elemToNode(ielem,5) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*j + i + 1;
297 elemToNode(ielem,6) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*(j + 1) + i + 1;
298 elemToNode(ielem,7) = (NY + 1)*(NX + 1)*(k + 1) + (NX + 1)*(j + 1) + i;
299 ielem++;
300 }
301 }
302 }
303
304 // Get edge connectivity
305 FieldContainer<int> edgeToNode(numEdges, numNodesPerEdge);
306 FieldContainer<int> elemToEdge(numElems, numEdgesPerElem);
307 int iedge = 0;
308 inode = 0;
309 for (int k=0; k<NZ+1; k++) {
310 for (int j=0; j<NY+1; j++) {
311 for (int i=0; i<NX+1; i++) {
312 if (i < NX){
313 edgeToNode(iedge,0) = inode;
314 edgeToNode(iedge,1) = inode + 1;
315 if (j < NY && k < NZ){
316 ielem=i+j*NX+k*NX*NY;
317 elemToEdge(ielem,0) = iedge;
318 if (j > 0)
319 elemToEdge(ielem-NX,2) = iedge;
320 if (k > 0)
321 elemToEdge(ielem-NX*NY,4) = iedge;
322 if (j > 0 && k > 0)
323 elemToEdge(ielem-NX*NY-NX,6) = iedge;
324 }
325 else if (j == NY && k == NZ){
326 ielem=i+(NY-1)*NX+(NZ-1)*NX*NY;
327 elemToEdge(ielem,6) = iedge;
328 }
329 else if (k == NZ && j < NY){
330 ielem=i+j*NX+(NZ-1)*NX*NY;
331 elemToEdge(ielem,4) = iedge;
332 if (j > 0)
333 elemToEdge(ielem-NX,6) = iedge;
334 }
335 else if (k < NZ && j == NY){
336 ielem=i+(NY-1)*NX+k*NX*NY;
337 elemToEdge(ielem,2) = iedge;
338 if (k > 0)
339 elemToEdge(ielem-NX*NY,6) = iedge;
340 }
341 iedge++;
342 }
343 if (j < NY){
344 edgeToNode(iedge,0) = inode;
345 edgeToNode(iedge,1) = inode + NX+1;
346 if (i < NX && k < NZ){
347 ielem=i+j*NX+k*NX*NY;
348 elemToEdge(ielem,3) = iedge;
349 if (i > 0)
350 elemToEdge(ielem-1,1) = iedge;
351 if (k > 0)
352 elemToEdge(ielem-NX*NY,7) = iedge;
353 if (i > 0 && k > 0)
354 elemToEdge(ielem-NX*NY-1,5) = iedge;
355 }
356 else if (i == NX && k == NZ){
357 ielem=NX-1+j*NX+(NZ-1)*NX*NY;
358 elemToEdge(ielem,5) = iedge;
359 }
360 else if (k == NZ && i < NX){
361 ielem=i+j*NX+(NZ-1)*NX*NY;
362 elemToEdge(ielem,7) = iedge;
363 if (i > 0)
364 elemToEdge(ielem-1,5) = iedge;
365 }
366 else if (k < NZ && i == NX){
367 ielem=NX-1+j*NX+k*NX*NY;
368 elemToEdge(ielem,1) = iedge;
369 if (k > 0)
370 elemToEdge(ielem-NX*NY,5) = iedge;
371 }
372 iedge++;
373 }
374 if (k < NZ){
375 edgeToNode(iedge,0) = inode;
376 edgeToNode(iedge,1) = inode + (NX+1)*(NY+1);
377 if (i < NX && j < NY){
378 ielem=i+j*NX+k*NX*NY;
379 elemToEdge(ielem,8) = iedge;
380 if (i > 0)
381 elemToEdge(ielem-1,9) = iedge;
382 if (j > 0)
383 elemToEdge(ielem-NX,11) = iedge;
384 if (i > 0 && j > 0)
385 elemToEdge(ielem-NX-1,10) = iedge;
386 }
387 else if (i == NX && j == NY){
388 ielem=NX-1+(NY-1)*NX+k*NX*NY;
389 elemToEdge(ielem,10) = iedge;
390 }
391 else if (j == NY && i < NX){
392 ielem=i+(NY-1)*NX+k*NX*NY;
393 elemToEdge(ielem,11) = iedge;
394 if (i > 0)
395 elemToEdge(ielem-1,10) = iedge;
396 }
397 else if (j < NY && i == NX){
398 ielem=NX-1+j*NX+k*NX*NY;
399 elemToEdge(ielem,9) = iedge;
400 if (j > 0)
401 elemToEdge(ielem-NX,10) = iedge;
402 }
403 iedge++;
404 }
405 inode++;
406 }
407 }
408 }
409
410 // Find boundary edges
411 FieldContainer<int> edgeOnBoundary(numEdges);
412 for (int i=0; i<numEdges; i++){
413 if (nodeOnBoundary(edgeToNode(i,0)) && nodeOnBoundary(edgeToNode(i,1))){
414 edgeOnBoundary(i)=1;
415 }
416 }
417
418 // Get face connectivity
419 FieldContainer<int> faceToNode(numFaces, numNodesPerFace);
420 FieldContainer<int> elemToFace(numElems, numFacesPerElem);
421 FieldContainer<int> faceToEdge(numFaces, numEdgesPerFace);
422 int iface = 0;
423 inode = 0;
424 for (int k=0; k<NZ+1; k++) {
425 for (int j=0; j<NY+1; j++) {
426 for (int i=0; i<NX+1; i++) {
427 if (i < NX && k < NZ) {
428 faceToNode(iface,0)=inode;
429 faceToNode(iface,1)=inode + 1;
430 faceToNode(iface,2)=inode + (NX+1)*(NY+1)+1;
431 faceToNode(iface,3)=inode + (NX+1)*(NY+1);
432 if (j < NY) {
433 ielem=i+j*NX+k*NX*NY;
434 faceToEdge(iface,0)=elemToEdge(ielem,0);
435 faceToEdge(iface,1)=elemToEdge(ielem,9);
436 faceToEdge(iface,2)=elemToEdge(ielem,4);
437 faceToEdge(iface,3)=elemToEdge(ielem,8);
438 elemToFace(ielem,0)=iface;
439 if (j > 0) {
440 elemToFace(ielem-NX,2)=iface;
441 }
442 }
443 else if (j == NY) {
444 ielem=i+(NY-1)*NX+k*NX*NY;
445 faceToEdge(iface,0)=elemToEdge(ielem,2);
446 faceToEdge(iface,1)=elemToEdge(ielem,10);
447 faceToEdge(iface,2)=elemToEdge(ielem,6);
448 faceToEdge(iface,3)=elemToEdge(ielem,11);
449 elemToFace(ielem,2)=iface;
450 }
451 iface++;
452 }
453 if (j < NY && k < NZ) {
454 faceToNode(iface,0)=inode;
455 faceToNode(iface,1)=inode + NX+1;
456 faceToNode(iface,2)=inode + (NX+1)*(NY+1) + NX+1;
457 faceToNode(iface,3)=inode + (NX+1)*(NY+1);
458 if (i < NX) {
459 ielem=i+j*NX+k*NX*NY;
460 faceToEdge(iface,0)=elemToEdge(ielem,3);
461 faceToEdge(iface,1)=elemToEdge(ielem,11);
462 faceToEdge(iface,2)=elemToEdge(ielem,7);
463 faceToEdge(iface,3)=elemToEdge(ielem,8);
464 elemToFace(ielem,3)=iface;
465 if (i > 0) {
466 elemToFace(ielem-1,1)=iface;
467 }
468 }
469 else if (i == NX) {
470 ielem=NX-1+j*NX+k*NX*NY;
471 faceToEdge(iface,0)=elemToEdge(ielem,1);
472 faceToEdge(iface,1)=elemToEdge(ielem,10);
473 faceToEdge(iface,2)=elemToEdge(ielem,5);
474 faceToEdge(iface,3)=elemToEdge(ielem,9);
475 elemToFace(ielem,1)=iface;
476 }
477 iface++;
478 }
479 if (i < NX && j < NY) {
480 faceToNode(iface,0)=inode;
481 faceToNode(iface,1)=inode + 1;
482 faceToNode(iface,2)=inode + NX+2;
483 faceToNode(iface,3)=inode + NX+1;
484 if (k < NZ) {
485 ielem=i+j*NX+k*NX*NY;
486 faceToEdge(iface,0)=elemToEdge(ielem,0);
487 faceToEdge(iface,1)=elemToEdge(ielem,1);
488 faceToEdge(iface,2)=elemToEdge(ielem,2);
489 faceToEdge(iface,3)=elemToEdge(ielem,3);
490 elemToFace(ielem,4)=iface;
491 if (k > 0) {
492 elemToFace(ielem-NX*NY,5)=iface;
493 }
494 }
495 else if (k == NZ) {
496 ielem=i+j*NX+(NZ-1)*NX*NY;
497 faceToEdge(iface,0)=elemToEdge(ielem,4);
498 faceToEdge(iface,1)=elemToEdge(ielem,5);
499 faceToEdge(iface,2)=elemToEdge(ielem,6);
500 faceToEdge(iface,3)=elemToEdge(ielem,7);
501 elemToFace(ielem,5)=iface;
502 }
503 iface++;
504 }
505 inode++;
506 }
507 }
508 }
509
510 // Find boundary faces
511 FieldContainer<int> faceOnBoundary(numFaces);
512 for (int i=0; i<numFaces; i++){
513 if (nodeOnBoundary(faceToNode(i,0)) && nodeOnBoundary(faceToNode(i,1))
514 && nodeOnBoundary(faceToNode(i,2)) && nodeOnBoundary(faceToNode(i,3))){
515 faceOnBoundary(i)=1;
516 }
517 }
518
519#define DUMP_DATA
520#ifdef DUMP_DATA
521 // Output connectivity
522 ofstream fe2nout("elem2node.dat");
523 ofstream fe2fout("elem2face.dat");
524 for (int k=0; k<NZ; k++) {
525 for (int j=0; j<NY; j++) {
526 for (int i=0; i<NX; i++) {
527 ielem = i + j * NX + k * NX * NY;
528 for (int m=0; m<numNodesPerElem; m++){
529 fe2nout << elemToNode(ielem,m) <<" ";
530 }
531 fe2nout <<"\n";
532 for (int n=0; n<numFacesPerElem; n++) {
533 fe2fout << elemToFace(ielem,n) << " ";
534 }
535 fe2fout << "\n";
536 }
537 }
538 }
539 fe2nout.close();
540 fe2fout.close();
541#endif
542
543#ifdef DUMP_DATA_EXTRA
544 ofstream ff2nout("face2node.dat");
545 ofstream ff2eout("face2edge.dat");
546 for (int i=0; i<numFaces; i++) {
547 for (int j=0; j<numNodesPerFace; j++) {
548 ff2nout << faceToNode(i,j) << " ";
549 }
550 for (int k=0; k<numEdgesPerFace; k++) {
551 ff2eout << faceToEdge(i,k) << " ";
552 }
553 ff2nout << "\n";
554 ff2eout << "\n";
555 }
556 ff2nout.close();
557 ff2eout.close();
558
559 ofstream fBnodeout("nodeOnBndy.dat");
560 ofstream fBfaceout("faceOnBndy.dat");
561 for (int i=0; i<numNodes; i++) {
562 fBnodeout << nodeOnBoundary(i) <<"\n";
563 }
564 for (int i=0; i<numFaces; i++) {
565 fBfaceout << faceOnBoundary(i) <<"\n";
566 }
567 fBnodeout.close();
568 fBfaceout.close();
569#endif
570
571 // Set material properties using undeformed grid assuming each element has only one value of mu
572 FieldContainer<double> muVal(numElems);
573 for (int k=0; k<NZ; k++) {
574 for (int j=0; j<NY; j++) {
575 for (int i=0; i<NX; i++) {
576 ielem = i + j * NX + k * NX * NY;
577 double midElemX = nodeCoord(elemToNode(ielem,0),0) + hx/2.0;
578 double midElemY = nodeCoord(elemToNode(ielem,0),1) + hy/2.0;
579 double midElemZ = nodeCoord(elemToNode(ielem,0),2) + hz/2.0;
580 if ( (midElemX > mu1LeftX) && (midElemY > mu1LeftY) && (midElemZ > mu1LeftZ) &&
581 (midElemX <= mu1RightX) && (midElemY <= mu1RightY) && (midElemZ <= mu1RightZ) ){
582 muVal(ielem) = mu1;
583 }
584 else {
585 muVal(ielem) = mu2;
586 }
587 }
588 }
589 }
590
591 // Perturb mesh coordinates (only interior nodes)
592 if (randomMesh){
593 for (int k=1; k<NZ; k++) {
594 for (int j=1; j<NY; j++) {
595 for (int i=1; i<NX; i++) {
596 int inode = i + j * (NX + 1) + k * (NX + 1) * (NY + 1);
597 // random numbers between -1.0 and 1.0
598 double rx = 2.0 * (double)rand()/RAND_MAX - 1.0;
599 double ry = 2.0 * (double)rand()/RAND_MAX - 1.0;
600 double rz = 2.0 * (double)rand()/RAND_MAX - 1.0;
601 // limit variation to 1/4 edge length
602 nodeCoord(inode,0) = nodeCoord(inode,0) + 0.125 * hx * rx;
603 nodeCoord(inode,1) = nodeCoord(inode,1) + 0.125 * hy * ry;
604 nodeCoord(inode,2) = nodeCoord(inode,2) + 0.125 * hz * rz;
605 }
606 }
607 }
608 }
609
610#ifdef DUMP_DATA
611 // Print nodal coords
612 ofstream fcoordout("coords.dat");
613 for (int i=0; i<numNodes; i++) {
614 fcoordout << nodeCoord(i,0) <<" ";
615 fcoordout << nodeCoord(i,1) <<" ";
616 fcoordout << nodeCoord(i,2) <<"\n";
617 }
618 fcoordout.close();
619#endif
620
621
622// **************************** INCIDENCE MATRIX **************************************
623
624 // Edge to face incidence matrix
625 *outStream << "Building incidence matrix ... \n\n";
626
627 Epetra_SerialComm Comm;
628 Epetra_Map globalMapD(numFaces, 0, Comm);
629 Epetra_Map globalMapC(numEdges, 0, Comm);
630 Epetra_FECrsMatrix DCurl(Copy, globalMapD, globalMapC, 2);
631
632 double vals[4];
633 vals[0]=1.0; vals[1]=1.0; vals[2]=-1.0; vals[3]=-1.0;
634 for (int j=0; j<numFaces; j++){
635 int rowNum = j;
636 int colNum[4];
637 colNum[0] = faceToEdge(j,0);
638 colNum[1] = faceToEdge(j,1);
639 colNum[2] = faceToEdge(j,2);
640 colNum[3] = faceToEdge(j,3);
641 DCurl.InsertGlobalValues(1, &rowNum, 4, colNum, vals);
642 }
643
644
645// ************************************ CUBATURE **************************************
646
647 // Get numerical integration points and weights
648 *outStream << "Getting cubature ... \n\n";
649
651 int cubDegree = 2;
652 Teuchos::RCP<Cubature<double> > hexCub = cubFactory.create(hex_8, cubDegree);
653
654 int cubDim = hexCub->getDimension();
655 int numCubPoints = hexCub->getNumPoints();
656
657 FieldContainer<double> cubPoints(numCubPoints, cubDim);
658 FieldContainer<double> cubWeights(numCubPoints);
659
660 hexCub->getCubature(cubPoints, cubWeights);
661
662
663 // Get numerical integration points and weights for hexahedron face
664 // (needed for rhs boundary term)
665
666 // Define topology of the face parametrization domain as [-1,1]x[-1,1]
667 CellTopology paramQuadFace(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
668
669 // Define cubature
670 DefaultCubatureFactory<double> cubFactoryFace;
671 Teuchos::RCP<Cubature<double> > hexFaceCubature = cubFactoryFace.create(paramQuadFace, 3);
672 int cubFaceDim = hexFaceCubature -> getDimension();
673 int numFacePoints = hexFaceCubature -> getNumPoints();
674
675 // Define storage for cubature points and weights on [-1,1]x[-1,1]
676 FieldContainer<double> paramGaussWeights(numFacePoints);
677 FieldContainer<double> paramGaussPoints(numFacePoints,cubFaceDim);
678
679 // Define storage for cubature points on workset faces
680 hexFaceCubature -> getCubature(paramGaussPoints, paramGaussWeights);
681
682
683
684// ************************************** BASIS ***************************************
685
686 // Define basis
687 *outStream << "Getting basis ... \n\n";
690
691 int numFieldsC = hexHCurlBasis.getCardinality();
692 int numFieldsD = hexHDivBasis.getCardinality();
693
694 // Evaluate basis at cubature points
695 FieldContainer<double> hexCVals(numFieldsC, numCubPoints, spaceDim);
696 FieldContainer<double> hexDVals(numFieldsD, numCubPoints, spaceDim);
697 FieldContainer<double> hexDivs(numFieldsD, numCubPoints);
698 FieldContainer<double> worksetDVals(numFieldsD, numFacePoints, spaceDim);
699
700 hexHCurlBasis.getValues(hexCVals, cubPoints, OPERATOR_VALUE);
701 hexHDivBasis.getValues(hexDVals, cubPoints, OPERATOR_VALUE);
702 hexHDivBasis.getValues(hexDivs, cubPoints, OPERATOR_DIV);
703
704
705// ******** LOOP OVER ELEMENTS TO CREATE LOCAL MASS and STIFFNESS MATRICES *************
706
707
708 *outStream << "Building mass and stiffness matrices ... \n\n";
709
710 // Settings and data structures for mass and stiffness matrices
712 typedef FunctionSpaceTools fst;
713 int numCells = 1;
714
715 // Containers for nodes, edge and face signs
716 FieldContainer<double> hexNodes(numCells, numNodesPerElem, spaceDim);
717 FieldContainer<double> hexEdgeSigns(numCells, numFieldsC);
718 FieldContainer<double> hexFaceSigns(numCells, numFieldsD);
719 // Containers for Jacobian
720 FieldContainer<double> hexJacobian(numCells, numCubPoints, spaceDim, spaceDim);
721 FieldContainer<double> hexJacobInv(numCells, numCubPoints, spaceDim, spaceDim);
722 FieldContainer<double> hexJacobDet(numCells, numCubPoints);
723 // Containers for element HCURL mass matrix
724 FieldContainer<double> massMatrixC(numCells, numFieldsC, numFieldsC);
725 FieldContainer<double> weightedMeasure(numCells, numCubPoints);
726 FieldContainer<double> weightedMeasureMuInv(numCells, numCubPoints);
727 FieldContainer<double> hexCValsTransformed(numCells, numFieldsC, numCubPoints, spaceDim);
728 FieldContainer<double> hexCValsTransformedWeighted(numCells, numFieldsC, numCubPoints, spaceDim);
729 // Containers for element HDIV mass matrix
730 FieldContainer<double> massMatrixD(numCells, numFieldsD, numFieldsD);
731 FieldContainer<double> weightedMeasureMu(numCells, numCubPoints);
732 FieldContainer<double> hexDValsTransformed(numCells, numFieldsD, numCubPoints, spaceDim);
733 FieldContainer<double> hexDValsTransformedWeighted(numCells, numFieldsD, numCubPoints, spaceDim);
734 // Containers for element HDIV stiffness matrix
735 FieldContainer<double> stiffMatrixD(numCells, numFieldsD, numFieldsD);
736 FieldContainer<double> hexDivsTransformed(numCells, numFieldsD, numCubPoints);
737 FieldContainer<double> hexDivsTransformedWeighted(numCells, numFieldsD, numCubPoints);
738 // Containers for right hand side vectors
739 FieldContainer<double> rhsDatag(numCells, numCubPoints, cubDim);
740 FieldContainer<double> rhsDatah(numCells, numCubPoints);
741 FieldContainer<double> gD(numCells, numFieldsD);
742 FieldContainer<double> hD(numCells, numFieldsD);
743 FieldContainer<double> gDBoundary(numCells, numFieldsD);
744 FieldContainer<double> refGaussPoints(numFacePoints,spaceDim);
745 FieldContainer<double> worksetGaussPoints(numCells,numFacePoints,spaceDim);
746 FieldContainer<double> worksetJacobians(numCells, numFacePoints, spaceDim, spaceDim);
747 FieldContainer<double> worksetJacobDet(numCells, numFacePoints);
748 FieldContainer<double> worksetFaceN(numCells, numFacePoints, spaceDim);
749 FieldContainer<double> worksetVFieldVals(numCells, numFacePoints, spaceDim);
750 FieldContainer<double> worksetDValsTransformed(numCells, numFieldsD, numFacePoints, spaceDim);
751 FieldContainer<double> curluFace(numCells, numFacePoints, spaceDim);
752 FieldContainer<double> worksetDataCrossField(numCells, numFieldsD, numFacePoints, spaceDim);
753 // Container for cubature points in physical space
754 FieldContainer<double> physCubPoints(numCells,numCubPoints, cubDim);
755
756
757 // Global arrays in Epetra format
758 Epetra_FECrsMatrix MassC(Copy, globalMapC, numFieldsC);
759 Epetra_FECrsMatrix MassD(Copy, globalMapD, numFieldsD);
760 Epetra_FECrsMatrix StiffD(Copy, globalMapD, numFieldsD);
761 Epetra_FEVector rhsD(globalMapD);
762
763#ifdef DUMP_DATA
764 ofstream fSignsout("faceSigns.dat");
765#endif
766
767 // *** Element loop ***
768 for (int k=0; k<numElems; k++) {
769
770 // Physical cell coordinates
771 for (int i=0; i<numNodesPerElem; i++) {
772 hexNodes(0,i,0) = nodeCoord(elemToNode(k,i),0);
773 hexNodes(0,i,1) = nodeCoord(elemToNode(k,i),1);
774 hexNodes(0,i,2) = nodeCoord(elemToNode(k,i),2);
775 }
776
777 // Face signs
778 for (int j=0; j<numFacesPerElem; j++) {
779 hexFaceSigns(0,j) = -1.0;
780 for (int i=0; i<numNodesPerFace; i++) {
781 int indf=i+1;
782 if (indf > numNodesPerFace) indf=0;
783 if (elemToNode(k,refFaceToNode(j,0))==faceToNode(elemToFace(k,j),i) &&
784 elemToNode(k,refFaceToNode(j,1))==faceToNode(elemToFace(k,j),indf))
785 hexFaceSigns(0,j) = 1.0;
786 }
787#ifdef DUMP_DATA
788 fSignsout << hexFaceSigns(0,j) << " ";
789#endif
790 }
791#ifdef DUMP_DATA
792 fSignsout << "\n";
793#endif
794
795 // Edge signs
796 for (int j=0; j<numEdgesPerElem; j++) {
797 if (elemToNode(k,refEdgeToNode(j,0))==edgeToNode(elemToEdge(k,j),0) &&
798 elemToNode(k,refEdgeToNode(j,1))==edgeToNode(elemToEdge(k,j),1))
799 hexEdgeSigns(0,j) = 1.0;
800 else
801 hexEdgeSigns(0,j) = -1.0;
802 }
803
804 // Compute cell Jacobians, their inverses and their determinants
805 CellTools::setJacobian(hexJacobian, cubPoints, hexNodes, hex_8);
806 CellTools::setJacobianInv(hexJacobInv, hexJacobian );
807 CellTools::setJacobianDet(hexJacobDet, hexJacobian );
808
809
810// ************************** Compute element HCurl mass matrices *******************************
811
812 // transform to physical coordinates
813 fst::HCURLtransformVALUE<double>(hexCValsTransformed, hexJacobInv,
814 hexCVals);
815 // compute weighted measure
816 fst::computeCellMeasure<double>(weightedMeasure, hexJacobDet, cubWeights);
817
818 // multiply by weighted measure
819 fst::multiplyMeasure<double>(hexCValsTransformedWeighted,
820 weightedMeasure, hexCValsTransformed);
821
822 // integrate to compute element mass matrix
823 fst::integrate<double>(massMatrixC,
824 hexCValsTransformed, hexCValsTransformedWeighted,
825 COMP_CPP);
826 // apply edge signs
827 fst::applyLeftFieldSigns<double>(massMatrixC, hexEdgeSigns);
828 fst::applyRightFieldSigns<double>(massMatrixC, hexEdgeSigns);
829
830 // assemble into global matrix
831 for (int row = 0; row < numFieldsC; row++){
832 for (int col = 0; col < numFieldsC; col++){
833 int rowIndex = elemToEdge(k,row);
834 int colIndex = elemToEdge(k,col);
835 double val = massMatrixC(0,row,col);
836 MassC.InsertGlobalValues(1, &rowIndex, 1, &colIndex, &val);
837 }
838 }
839
840// ************************** Compute element HDiv mass matrices *******************************
841
842 // transform to physical coordinates
843 fst::HDIVtransformVALUE<double>(hexDValsTransformed, hexJacobian, hexJacobDet,
844 hexDVals);
845
846 // multiply by weighted measure
847 fst::multiplyMeasure<double>(hexDValsTransformedWeighted,
848 weightedMeasure, hexDValsTransformed);
849
850 // integrate to compute element mass matrix
851 fst::integrate<double>(massMatrixD,
852 hexDValsTransformed, hexDValsTransformedWeighted,
853 COMP_CPP);
854 // apply face signs
855 fst::applyLeftFieldSigns<double>(massMatrixD, hexFaceSigns);
856 fst::applyRightFieldSigns<double>(massMatrixD, hexFaceSigns);
857
858 // assemble into global matrix
859 for (int row = 0; row < numFieldsD; row++){
860 for (int col = 0; col < numFieldsD; col++){
861 int rowIndex = elemToFace(k,row);
862 int colIndex = elemToFace(k,col);
863 double val = massMatrixD(0,row,col);
864 MassD.InsertGlobalValues(1, &rowIndex, 1, &colIndex, &val);
865 }
866 }
867
868// ************************ Compute element HDiv stiffness matrices *****************************
869
870 // transform to physical coordinates
871 fst::HDIVtransformDIV<double>(hexDivsTransformed, hexJacobDet,
872 hexDivs);
873
874 // multiply by weighted measure
875 fst::multiplyMeasure<double>(hexDivsTransformedWeighted,
876 weightedMeasure, hexDivsTransformed);
877
878 // integrate to compute element stiffness matrix
879 fst::integrate<double>(stiffMatrixD,
880 hexDivsTransformed, hexDivsTransformedWeighted,
881 COMP_CPP);
882
883 // apply face signs
884 fst::applyLeftFieldSigns<double>(stiffMatrixD, hexFaceSigns);
885 fst::applyRightFieldSigns<double>(stiffMatrixD, hexFaceSigns);
886
887 // assemble into global matrix
888 for (int row = 0; row < numFieldsD; row++){
889 for (int col = 0; col < numFieldsD; col++){
890 int rowIndex = elemToFace(k,row);
891 int colIndex = elemToFace(k,col);
892 double val = stiffMatrixD(0,row,col);
893 StiffD.InsertGlobalValues(1, &rowIndex, 1, &colIndex, &val);
894 }
895 }
896
897// ******************************* Build right hand side ************************************
898
899 // transform integration points to physical points
900 CellTools::mapToPhysicalFrame(physCubPoints, cubPoints, hexNodes, hex_8);
901
902 // evaluate right hand side functions at physical points
903 for (int nPt = 0; nPt < numCubPoints; nPt++){
904
905 double x = physCubPoints(0,nPt,0);
906 double y = physCubPoints(0,nPt,1);
907 double z = physCubPoints(0,nPt,2);
908 double du1, du2, du3;
909
910 evalCurlCurlu(du1, du2, du3, x, y, z);
911 rhsDatag(0,nPt,0) = du1;
912 rhsDatag(0,nPt,1) = du2;
913 rhsDatag(0,nPt,2) = du3;
914
915 rhsDatah(0,nPt) = evalDivu(x, y, z);
916 }
917
918 // integrate (curl g, w) term
919 fst::integrate<double>(gD, rhsDatag, hexDValsTransformedWeighted,
920 COMP_CPP);
921
922 // integrate (h,div w) term
923 fst::integrate<double>(hD, rhsDatah, hexDivsTransformedWeighted,
924 COMP_CPP);
925
926 // apply signs
927 fst::applyFieldSigns<double>(gD, hexFaceSigns);
928 fst::applyFieldSigns<double>(hD, hexFaceSigns);
929
930 // calculate boundary term, (g x w, n)_{\Gamma}
931 for (int i = 0; i < numFacesPerElem; i++){
932 if (faceOnBoundary(elemToFace(k,i))){
933
934 // map Gauss points on quad to reference face: paramGaussPoints -> refGaussPoints
936 paramGaussPoints,
937 2, i, hex_8);
938
939 // get basis values at points on reference cell
940 hexHDivBasis.getValues(worksetDVals, refGaussPoints, OPERATOR_VALUE);
941
942 // compute Jacobians at Gauss pts. on reference face for all parent cells
943 CellTools::setJacobian(worksetJacobians, refGaussPoints,
944 hexNodes, hex_8);
945 CellTools::setJacobianDet(worksetJacobDet, worksetJacobians);
946
947 // transform to physical coordinates
948 fst::HDIVtransformVALUE<double>(worksetDValsTransformed, worksetJacobians,
949 worksetJacobDet, worksetDVals);
950
951 // map Gauss points on quad from ref. face to face workset: refGaussPoints -> worksetGaussPoints
952 CellTools::mapToPhysicalFrame(worksetGaussPoints,
953 refGaussPoints,
954 hexNodes, hex_8);
955
956 // compute face normals
958 worksetJacobians,
959 i, hex_8);
960
961 // evaluate curl u at face points
962 for(int nPt = 0; nPt < numFacePoints; nPt++){
963
964 double x = worksetGaussPoints(0, nPt, 0);
965 double y = worksetGaussPoints(0, nPt, 1);
966 double z = worksetGaussPoints(0, nPt, 2);
967
968 evalCurlu(curluFace(0,nPt,0), curluFace(0,nPt,1), curluFace(0,nPt,2), x, y, z);
969 }
970
971 // compute the cross product of curluFace with basis and multiply by weights
972 for (int nF = 0; nF < numFieldsD; nF++){
973 for(int nPt = 0; nPt < numFacePoints; nPt++){
974 worksetDataCrossField(0,nF,nPt,0) = (curluFace(0,nPt,1)*worksetDValsTransformed(0,nF,nPt,2)
975 - curluFace(0,nPt,2)*worksetDValsTransformed(0,nF,nPt,1))
976 * paramGaussWeights(nPt);
977 worksetDataCrossField(0,nF,nPt,1) = (curluFace(0,nPt,2)*worksetDValsTransformed(0,nF,nPt,0)
978 - curluFace(0,nPt,0)*worksetDValsTransformed(0,nF,nPt,2))
979 * paramGaussWeights(nPt);
980 worksetDataCrossField(0,nF,nPt,2) = (curluFace(0,nPt,0)*worksetDValsTransformed(0,nF,nPt,1)
981 - curluFace(0,nPt,1)*worksetDValsTransformed(0,nF,nPt,0))
982 *paramGaussWeights(nPt);
983 } //nPt
984 } //nF
985
986 // Integrate
987 fst::integrate<double>(gDBoundary, worksetFaceN, worksetDataCrossField,
988 COMP_CPP);
989
990 // apply signs
991 fst::applyFieldSigns<double>(gDBoundary, hexFaceSigns);
992
993 // add into gD term
994 for (int nF = 0; nF < numFieldsD; nF++){
995 gD(0,nF) = gD(0,nF) - gDBoundary(0,nF);
996 }
997
998 } // if faceOnBoundary
999 } // numFaces
1000
1001
1002 // assemble into global vector
1003 for (int row = 0; row < numFieldsD; row++){
1004 int rowIndex = elemToFace(k,row);
1005 double val = hD(0,row)+gD(0,row);
1006 rhsD.SumIntoGlobalValues(1, &rowIndex, &val);
1007 }
1008
1009
1010 } // *** end element loop ***
1011
1012 // Assemble over multiple processors, if necessary
1013 DCurl.GlobalAssemble(); DCurl.FillComplete(MassC.RowMap(),MassD.RowMap());
1014 MassC.GlobalAssemble(); MassC.FillComplete();
1015 MassD.GlobalAssemble(); MassD.FillComplete();
1016 StiffD.GlobalAssemble(); StiffD.FillComplete();
1017 rhsD.GlobalAssemble();
1018
1019#ifdef DUMP_DATA
1020 fSignsout.close();
1021#endif
1022
1023 // Build the inverse diagonal for MassC
1024 Epetra_CrsMatrix MassCinv(Copy,MassC.RowMap(),MassC.RowMap(),1);
1025 Epetra_Vector DiagC(MassC.RowMap());
1026
1027 DiagC.PutScalar(1.0);
1028 MassC.Multiply(false,DiagC,DiagC);
1029 for(int i=0; i<DiagC.MyLength(); i++) {
1030 DiagC[i]=1.0/DiagC[i];
1031 }
1032 for(int i=0; i<DiagC.MyLength(); i++) {
1033 int CID=MassC.GCID(i);
1034 MassCinv.InsertGlobalValues(MassC.GRID(i),1,&(DiagC[i]),&CID);
1035 }
1036 MassCinv.FillComplete();
1037
1038 // Set value to zero on diagonal that cooresponds to boundary edge
1039 for(int i=0;i<numEdges;i++) {
1040 if (edgeOnBoundary(i)){
1041 double val=0.0;
1042 MassCinv.ReplaceGlobalValues(i,1,&val,&i);
1043 }
1044 }
1045
1046 int numEntries;
1047 double *values;
1048 int *cols;
1049
1050 // Adjust matrices and rhs due to boundary conditions
1051 for (int row = 0; row<numFaces; row++){
1052 MassD.ExtractMyRowView(row,numEntries,values,cols);
1053 for (int i=0; i<numEntries; i++){
1054 if (faceOnBoundary(cols[i])) {
1055 values[i]=0;
1056 }
1057 }
1058 StiffD.ExtractMyRowView(row,numEntries,values,cols);
1059 for (int i=0; i<numEntries; i++){
1060 if (faceOnBoundary(cols[i])) {
1061 values[i]=0;
1062 }
1063 }
1064 }
1065 for (int row = 0; row<numFaces; row++){
1066 if (faceOnBoundary(row)) {
1067 int rowindex = row;
1068 StiffD.ExtractMyRowView(row,numEntries,values,cols);
1069 for (int i=0; i<numEntries; i++){
1070 values[i]=0;
1071 }
1072 MassD.ExtractMyRowView(row,numEntries,values,cols);
1073 for (int i=0; i<numEntries; i++){
1074 values[i]=0;
1075 }
1076 rhsD[0][row]=0;
1077 double val = 1.0;
1078 StiffD.ReplaceGlobalValues(1, &rowindex, 1, &rowindex, &val);
1079 }
1080 }
1081
1082#ifdef DUMP_DATA
1083 // Dump matrices to disk
1084 EpetraExt::RowMatrixToMatlabFile("mag_m1inv_matrix.dat",MassCinv);
1085 EpetraExt::RowMatrixToMatlabFile("mag_m2_matrix.dat",MassD);
1086 EpetraExt::RowMatrixToMatlabFile("mag_k2_matrix.dat",StiffD);
1087 EpetraExt::RowMatrixToMatlabFile("mag_t1_matrix.dat",DCurl);
1088 EpetraExt::MultiVectorToMatrixMarketFile("mag_rhs2_vector.dat",rhsD,0,0,false);
1089#endif
1090
1091 std::cout << "End Result: TEST PASSED\n";
1092
1093 // reset format state of std::cout
1094 std::cout.copyfmt(oldFormatState);
1095
1096 return 0;
1097}
1098// Calculates value of exact solution u
1099 int evalu(double & uExact0, double & uExact1, double & uExact2, double & x, double & y, double & z)
1100 {
1101
1102 // function 1
1103 uExact0 = exp(y+z)*(x+1.0)*(x-1.0);
1104 uExact1 = exp(x+z)*(y+1.0)*(y-1.0);
1105 uExact2 = exp(x+y)*(z+1.0)*(z-1.0);
1106
1107 /*
1108 // function 2
1109 uExact0 = cos(M_PI*y)*cos(M_PI*z)*(x+1.0)*(x-1.0);
1110 uExact1 = cos(M_PI*x)*cos(M_PI*z)*(y+1.0)*(y-1.0);
1111 uExact2 = cos(M_PI*x)*cos(M_PI*y)*(z+1.0)*(z-1.0);
1112
1113 */
1114 /*
1115 // function 3
1116 uExact0 = x*x-1.0;
1117 uExact1 = y*y-1.0;
1118 uExact2 = z*z-1.0;
1119
1120 // function 4
1121 uExact0 = sin(M_PI*x);
1122 uExact1 = sin(M_PI*y);
1123 uExact2 = sin(M_PI*z);
1124 */
1125
1126 return 0;
1127 }
1128// Calculates divergence of exact solution u
1129 double evalDivu(double & x, double & y, double & z)
1130 {
1131
1132 // function 1
1133 double divu = 2.0*x*exp(y+z)+2.0*y*exp(x+z)+2.0*z*exp(x+y);
1134
1135 // function 2
1136 // double divu = 2.0*x*cos(M_PI*y)*cos(M_PI*z) + 2.0*y*cos(M_PI*x)*cos(M_PI*z)
1137 // + 2.0*z*cos(M_PI*x)*cos(M_PI*y);
1138
1139 // function 3
1140 // double divu = 2.0*(x + y + z);
1141
1142 // function 4
1143 // double divu = M_PI*(cos(M_PI*x)+cos(M_PI*y)+cos(M_PI*z));
1144
1145 return divu;
1146 }
1147// Calculates curl of exact solution u
1148 int evalCurlu(double & curlu0, double & curlu1, double & curlu2, double & x, double & y, double & z)
1149 {
1150
1151 // function 1
1152 double duxdy = exp(y+z)*(x+1.0)*(x-1.0);
1153 double duxdz = exp(y+z)*(x+1.0)*(x-1.0);
1154 double duydx = exp(x+z)*(y+1.0)*(y-1.0);
1155 double duydz = exp(x+z)*(y+1.0)*(y-1.0);
1156 double duzdx = exp(x+y)*(z+1.0)*(z-1.0);
1157 double duzdy = exp(x+y)*(z+1.0)*(z-1.0);
1158
1159 /*
1160
1161 // function 2
1162 double duxdy = -M_PI*sin(M_PI*y)*cos(M_PI*z)*(x+1.0)*(x-1.0);
1163 double duxdz = -M_PI*sin(M_PI*z)*cos(M_PI*y)*(x+1.0)*(x-1.0);
1164 double duydx = -M_PI*sin(M_PI*x)*cos(M_PI*z)*(y+1.0)*(y-1.0);
1165 double duydz = -M_PI*sin(M_PI*z)*cos(M_PI*x)*(y+1.0)*(y-1.0);
1166 double duzdx = -M_PI*sin(M_PI*x)*cos(M_PI*y)*(z+1.0)*(z-1.0);
1167 double duzdy = -M_PI*sin(M_PI*y)*cos(M_PI*x)*(z+1.0)*(z-1.0);
1168 */
1169
1170 curlu0 = duzdy - duydz;
1171 curlu1 = duxdz - duzdx;
1172 curlu2 = duydx - duxdy;
1173
1174 /*
1175 // function 3 and 4
1176 curlu0 = 0;
1177 curlu1 = 0;
1178 curlu2 = 0;
1179 */
1180
1181 return 0;
1182 }
1183
1184// Calculates curl of the curl of exact solution u
1185 int evalCurlCurlu(double & curlCurlu0, double & curlCurlu1, double & curlCurlu2, double & x, double & y, double & z)
1186{
1187
1188 // function 1
1189 double dcurlu0dy = exp(x+y)*(z+1.0)*(z-1.0) - 2.0*y*exp(x+z);
1190 double dcurlu0dz = 2.0*z*exp(x+y) - exp(x+z)*(y+1.0)*(y-1.0);
1191 double dcurlu1dx = 2.0*x*exp(y+z) - exp(x+y)*(z+1.0)*(z-1.0);
1192 double dcurlu1dz = exp(y+z)*(x+1.0)*(x-1.0) - 2.0*z*exp(x+y);
1193 double dcurlu2dx = exp(x+z)*(y+1.0)*(y-1.0) - 2.0*x*exp(y+z);
1194 double dcurlu2dy = 2.0*y*exp(x+z) - exp(y+z)*(x+1.0)*(x-1.0);
1195
1196 /*
1197 // function 2
1198 double dcurlu0dy = -M_PI*M_PI*cos(M_PI*y)*cos(M_PI*x)*(z+1.0)*(z-1.0)
1199 + 2.0*y*M_PI*sin(M_PI*z)*cos(M_PI*x);
1200 double dcurlu0dz = -2.0*z*M_PI*sin(M_PI*y)*cos(M_PI*x)
1201 + M_PI*M_PI*cos(M_PI*z)*cos(M_PI*x)*(y+1.0)*(y-1.0);
1202 double dcurlu1dx = -2.0*x*M_PI*sin(M_PI*z)*cos(M_PI*y)
1203 + M_PI*M_PI*cos(M_PI*x)*cos(M_PI*y)*(z+1.0)*(z-1.0);
1204 double dcurlu1dz = -M_PI*M_PI*cos(M_PI*z)*cos(M_PI*y)*(x+1.0)*(x-1.0)
1205 + 2.0*z*M_PI*sin(M_PI*x)*cos(M_PI*y);
1206 double dcurlu2dx = -M_PI*M_PI*cos(M_PI*x)*cos(M_PI*z)*(y+1.0)*(y-1.0)
1207 + 2.0*x*M_PI*sin(M_PI*y)*cos(M_PI*z);
1208 double dcurlu2dy = -2.0*y*M_PI*sin(M_PI*x)*cos(M_PI*z)
1209 + M_PI*M_PI*cos(M_PI*y)*cos(M_PI*z)*(x+1.0)*(x-1.0);
1210 */
1211
1212 curlCurlu0 = dcurlu2dy - dcurlu1dz;
1213 curlCurlu1 = dcurlu0dz - dcurlu2dx;
1214 curlCurlu2 = dcurlu1dx - dcurlu0dy;
1215
1216 /*
1217 // function 3 and 4
1218 curlCurlu0 = 0.0;
1219 curlCurlu1 = 0.0;
1220 curlCurlu2 = 0.0;
1221 */
1222
1223 return 0;
1224}
1225
1226
Header file for utility class to provide array tools, such as tensor contractions,...
Header file for the Intrepid::CellTools class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for the Intrepid::HCURL_HEX_I1_FEM class.
Header file for the Intrepid::HDIV_HEX_I1_FEM class.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Intrepid utilities.
Implementation of the default H(curl)-compatible FEM basis of degree 1 on Hexahedron cell.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedron cell.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
virtual int getCardinality() const
Returns cardinality of the basis.
A stateless class for operations on cell data. Provides methods for:
static void getPhysicalFaceNormals(ArrayFaceNormal &faceNormals, const ArrayJac &worksetJacobians, const int &worksetFaceOrd, const shards::CellTopology &parentCell)
Computes non-normalized normal vectors to physical faces in a face workset ; (see Subcell worksets fo...
static void mapToReferenceSubcell(ArraySubcellPoint &refSubcellPoints, const ArrayParamPoint &paramPoints, const int subcellDim, const int subcellOrd, const shards::CellTopology &parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
static void mapToPhysicalFrame(ArrayPhysPoint &physPoints, const ArrayRefPoint &refPoints, const ArrayCell &cellWorkset, const shards::CellTopology &cellTopo, const int &whichCell=-1)
Computes F, the reference-to-physical frame map.
static void setJacobianDet(ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F.
static void setJacobianInv(ArrayJacInv &jacobianInv, const ArrayJac &jacobian)
Computes the inverse of the Jacobian matrix DF of the reference-to-physical frame map F.
A factory class that generates specific instances of cubatures.
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > &degree)
Factory method.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
Defines expert-level interfaces for the evaluation of functions and operators in physical space (supp...