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