Anasazi Version of the Day
Loading...
Searching...
No Matches
ModeLaplace3DQ2.cpp
1// @HEADER
2// *****************************************************************************
3// Anasazi: Block Eigensolvers Package
4//
5// Copyright 2004 NTESS and the Anasazi contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10//**************************************************************************
11//
12// NOTICE
13//
14// This software is a result of the research described in the report
15//
16// " A comparison of algorithms for modal analysis in the absence
17// of a sparse direct method", P. ArbenZ, R. Lehoucq, and U. Hetmaniuk,
18// Sandia National Laboratories, Technical report SAND2003-1028J.
19//
20// It is based on the Epetra, AztecOO, and ML packages defined in the Trilinos
21// framework ( http://software.sandia.gov/trilinos/ ).
22//
23// The distribution of this software follows also the rules defined in Trilinos.
24// This notice shall be marked on anY reproduction of this software, in whole or
25// in part.
26//
27// This program is distributed in the hope that it will be useful, but
28// WITHOUT ANY WARRANTY; without even the implied warranty of
29// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
30//
31// Code Authors: U. Hetmaniuk (ulhetma@sandia.gov), R. Lehoucq (rblehou@sandia.gov)
32//
33//**************************************************************************
34
35#include "ModeLaplace3DQ2.h"
36#include "Teuchos_Assert.hpp"
37
38
39const int ModeLaplace3DQ2::dofEle = 27;
40const int ModeLaplace3DQ2::maxConnect = 125;
41#ifndef M_PI
42const double ModeLaplace3DQ2::M_PI = 3.14159265358979323846;
43#endif
44
45
46ModeLaplace3DQ2::ModeLaplace3DQ2(const Epetra_Comm &_Comm, double _Lx, int _nX,
47 double _Ly, int _nY, double _Lz, int _nZ)
48 : myVerify(_Comm),
49 MyComm(_Comm),
50 mySort(),
51 Map(0),
52 K(0),
53 M(0),
54 Lx(_Lx),
55 nX(_nX),
56 Ly(_Ly),
57 nY(_nY),
58 Lz(_Lz),
59 nZ(_nZ),
60 x(0),
61 y(0),
62 z(0)
63 {
64
65 preProcess();
66
67}
68
69
70ModeLaplace3DQ2::~ModeLaplace3DQ2() {
71
72 if (Map)
73 delete Map;
74 Map = 0;
75
76 if (K)
77 delete K;
78 K = 0;
79
80 if (M)
81 delete M;
82 M = 0;
83
84 if (x)
85 delete[] x;
86 x = 0;
87
88 if (y)
89 delete[] y;
90 y = 0;
91
92 if (z)
93 delete[] z;
94 z = 0;
95
96}
97
98
99void ModeLaplace3DQ2::preProcess() {
100
101 // Create the distribution of equations across processors
102 makeMap();
103
104 // Count the number of elements touched by this processor
105 bool *isTouched = new bool[nX*nY*nZ];
106 int i, j, k;
107 for (k = 0; k < nZ; ++k)
108 for (j = 0; j < nY; ++j)
109 for (i=0; i<nX; ++i)
110 isTouched[i + j*nX + k*nX*nY] = false;
111
112 int numEle = countElements(isTouched);
113
114 // Create the mesh
115 int *elemTopo = new int[dofEle*numEle];
116 makeMyElementsTopology(elemTopo, isTouched);
117
118 delete[] isTouched;
119
120 // Get the number of nonZeros per row
121 int localSize = Map->NumMyElements();
122 int *numNz = new int[localSize];
123 int *connectivity = new int[localSize*maxConnect];
124 makeMyConnectivity(elemTopo, numEle, connectivity, numNz);
125
126 // Make the stiffness matrix
127 makeStiffness(elemTopo, numEle, connectivity, numNz);
128
129 // Assemble the mass matrix
130 makeMass(elemTopo, numEle, connectivity, numNz);
131
132 // Free some memory
133 delete[] elemTopo;
134 delete[] numNz;
135 delete[] connectivity;
136
137 // Get the geometrical coordinates of the managed nodes
138 double hx = Lx/nX;
139 double hy = Ly/nY;
140 double hz = Lz/nZ;
141
142 x = new double[localSize];
143 y = new double[localSize];
144 z = new double[localSize];
145
146 for (k=0; k<2*nZ-1; ++k) {
147 for (j=0; j<2*nY-1; ++j) {
148 for (i=0; i<2*nX-1; ++i) {
149 int node = i + j*(2*nX-1) + k*(2*nX-1)*(2*nY-1);
150 if (Map->LID(node) > -1) {
151 x[Map->LID(node)] = (i+1)*hx*0.5;
152 y[Map->LID(node)] = (j+1)*hy*0.5;
153 z[Map->LID(node)] = (k+1)*hz*0.5;
154 }
155 }
156 }
157 }
158
159}
160
161
162void ModeLaplace3DQ2::makeMap() {
163
164 int numProc = MyComm.NumProc();
165 int globalSize = (2*nX - 1)*(2*nY - 1)*(2*nZ - 1);
166 assert(globalSize > numProc);
167
168#ifdef _USE_CHACO
169 // Use the partitioner Chaco to distribute the unknowns
170 int *start = new int[globalSize+1];
171 memset(start, 0, (globalSize+1)*sizeof(int));
172
173 int i, j, k;
174 for (k=0; k<2*nZ-1; ++k) {
175 for (j=0; j<2*nY-1; ++j) {
176 for (i=0; i<2*nX-1; ++i) {
177 int node = i + j*(2*nX-1) + k*(2*nX-1)*(2*nY-1);
178 int connect = 1;
179 if (i%2 == 1) {
180 connect *= ((i == 1) || (i == 2*nX-3)) ? 4 : 5;
181 }
182 else {
183 connect *= ((i == 0) || (i == 2*nX-2)) ? 2 : 3;
184 }
185 if (j%2 == 1) {
186 connect *= ((j == 1) || (j == 2*nY-3)) ? 4 : 5;
187 }
188 else {
189 connect *= ((j == 0) || (j == 2*nY-2)) ? 2 : 3;
190 }
191 if (k%2 == 1) {
192 connect *= ((k == 1) || (k == 2*nZ-3)) ? 4 : 5;
193 }
194 else {
195 connect *= ((k == 0) || (k == 2*nZ-2)) ? 2 : 3;
196 }
197 // Don't count the node itself for Chaco
198 start[node+1] = connect - 1;
199 }
200 }
201 }
202
203 for (i = 0; i < globalSize; ++i)
204 start[i+1] += start[i];
205
206 int *adjacency = new int[start[globalSize]];
207 memset(adjacency, 0, start[globalSize]*sizeof(int));
208
209 int *elemTopo = new int[dofEle];
210 for (k=0; k<nZ; ++k) {
211 for (j=0; j<nY; ++j) {
212 for (i=0; i<nX; ++i) {
213 int middleMiddle = 2*i + 2*j*(2*nX - 1) + 2*k*(2*nX-1)*(2*nY-1);
214 elemTopo[26] = middleMiddle;
215 elemTopo[25] = (i == nX-1) ? -1 : middleMiddle + 1;
216 elemTopo[24] = (i == 0) ? -1 : middleMiddle - 1;
217 elemTopo[23] = (j == nY-1) ? -1 : middleMiddle + 2*nX - 1;
218 elemTopo[22] = (j == 0) ? -1 : middleMiddle - 2*nX + 1;
219 elemTopo[21] = (k == nZ-1) ? -1 : middleMiddle + (2*nX-1)*(2*nY-1) ;
220 elemTopo[20] = (k == 0) ? -1 : middleMiddle - (2*nX-1)*(2*nY-1);
221 elemTopo[19] = ((i == 0) || (j == nY-1)) ? -1 :
222 elemTopo[23] - 1;
223 elemTopo[18] = ((i == nX-1) || (j == nY-1)) ? -1 :
224 elemTopo[23] + 1;
225 elemTopo[17] = ((i == nX-1) || (j == 0)) ? -1 :
226 elemTopo[22] + 1;
227 elemTopo[16] = ((i == 0) || (j == 0)) ? -1 :
228 elemTopo[22] - 1;
229 elemTopo[15] = ((i == 0) || (k == nZ-1)) ? -1 :
230 elemTopo[21] - 1;
231 elemTopo[14] = ((j == nY-1) || (k == nZ-1)) ? -1 :
232 elemTopo[21] + 2*nX - 1;
233 elemTopo[13] = ((i == nX-1) || (k == nZ-1)) ? -1 :
234 elemTopo[21] + 1;
235 elemTopo[12] = ((j == 0) || (k == nZ-1)) ? -1 :
236 elemTopo[21] - 2*nX + 1;
237 elemTopo[11] = ((i == 0) || (k == 0)) ? -1 :
238 elemTopo[20] - 1;
239 elemTopo[10] = ((j == nY-1) || (k == 0)) ? -1 :
240 elemTopo[20] + 2*nX - 1;
241 elemTopo[ 9] = ((i == nX-1) || (k == 0)) ? -1 :
242 elemTopo[20] + 1;
243 elemTopo[ 8] = ((j == 0) || (k == 0)) ? -1 :
244 elemTopo[20] - 2*nX + 1;
245 elemTopo[ 7] = ((i == 0) || (j == nY-1) || (k == nZ-1)) ? -1 :
246 elemTopo[21] + 2*nX - 2;
247 elemTopo[ 6] = ((i == nX-1) || (j == nY-1) || (k == nZ-1)) ? -1 :
248 elemTopo[21] + 2*nX;
249 elemTopo[ 5] = ((i == nX-1) || (j == 0) || (k == nZ-1)) ? -1 :
250 elemTopo[21] - 2*nX + 2;
251 elemTopo[ 4] = ((i == 0) || (j == 0) || (k == nZ-1)) ? -1 :
252 elemTopo[21] - 2*nX;
253 elemTopo[ 3] = ((i == 0) || (j == nY-1) || (k == 0)) ? -1 :
254 elemTopo[20] + 2*nX - 2;
255 elemTopo[ 2] = ((i == nX-1) || (j == nY-1) || (k == 0)) ? -1 :
256 elemTopo[20] + 2*nX;
257 elemTopo[ 1] = ((i == nX-1) || (j == 0) || (k == 0)) ? -1 :
258 elemTopo[20] - 2*nX + 2;
259 elemTopo[ 0] = ((i == 0) || (j == 0) || (k == 0)) ? -1 :
260 elemTopo[20] - 2*nX;
261 for (int iD = 0; iD < dofEle; ++iD) {
262 if (elemTopo[iD] == -1)
263 continue;
264 for (int jD = 0; jD < dofEle; ++jD) {
265 int neighbor = elemTopo[jD];
266 // Don't count the node itself for Chaco
267 if ((neighbor == -1) || (iD == jD))
268 continue;
269 // Check if this neighbor is already stored
270 for (int l = start[elemTopo[iD]]; l < start[elemTopo[iD]+1]; ++l) {
271 // Note that Chaco uses a Fortran numbering
272 if (adjacency[l] == neighbor + 1)
273 break;
274 if (adjacency[l] == 0) {
275 // Store the node ID
276 // Note that Chaco uses a Fortran numbering
277 adjacency[l] = elemTopo[jD] + 1;
278 break;
279 }
280 } // for (int l = start[elemTopo[iD]]; l < start[elemTopo[iD]+1]; ++l)
281 } // for (int jD = 0; jD < dofEle; ++jD)
282 } // for (int iD = 0; iD < dofEle; ++iD)
283 }
284 }
285 }
286 delete[] elemTopo;
287
288 int nDir[3];
289 nDir[0] = numProc;
290 nDir[1] = 0;
291 nDir[2] = 0;
292 short int *partition = new short int[globalSize];
293 // Call Chaco to partition the matrix
294 interface(globalSize, start, adjacency, 0, 0, 0, 0, 0, 0, 0,
295 partition, 1, 1, nDir, 0, 1, 1, 1, 250, 1, 0.001, 7654321L);
296 // Define the Epetra_Map
297 int localSize = 0;
298 int myPid = MyComm.MyPID();
299 for (i = 0; i < globalSize; ++i) {
300 if (partition[i] == myPid)
301 localSize += 1;
302 }
303 int *myRows = new int[localSize];
304 localSize = 0;
305 for (i = 0; i < globalSize; ++i) {
306 if (partition[i] == myPid) {
307 myRows[localSize] = i;
308 localSize +=1 ;
309 }
310 }
311 delete[] partition;
312 delete[] adjacency;
313 delete[] start;
314 Map = new Epetra_Map(globalSize, localSize, myRows, 0, MyComm);
315 delete[] myRows;
316#else
317 // Create a uniform distribution of the unknowns across processors
318 Map = new Epetra_Map(globalSize, 0, MyComm);
319#endif
320
321}
322
323
324int ModeLaplace3DQ2::countElements(bool *isTouched) {
325
326 // This routine counts and flags the elements that contain the nodes
327 // on this processor.
328
329 int i, j, k;
330 int numEle = 0;
331
332 for (k=0; k<nZ; ++k) {
333 for (j=0; j<nY; ++j) {
334 for (i=0; i<nX; ++i) {
335 isTouched[i + j*nX + k*nX*nY] = false;
336 int middleMiddle = 2*i + 2*j*(2*nX - 1) + 2*k*(2*nX-1)*(2*nY-1);
337 int node;
338 node = middleMiddle;
339 if ((node > -1) && (Map->LID(node) > -1)) {
340 isTouched[i + j*nX + k*nX*nY] = true;
341 numEle += 1;
342 continue;
343 }
344 node = (i == nX-1) ? -1 : middleMiddle + 1;
345 if ((node > -1) && (Map->LID(node) > -1)) {
346 isTouched[i + j*nX + k*nX*nY] = true;
347 numEle += 1;
348 continue;
349 }
350 node = (i == 0) ? -1 : middleMiddle - 1;
351 if ((node > -1) && (Map->LID(node) > -1)) {
352 isTouched[i + j*nX + k*nX*nY] = true;
353 numEle += 1;
354 continue;
355 }
356 node = (j == nY-1) ? -1 : middleMiddle + 2*nX - 1;
357 if ((node > -1) && (Map->LID(node) > -1)) {
358 isTouched[i + j*nX + k*nX*nY] = true;
359 numEle += 1;
360 continue;
361 }
362 node = (j == 0) ? -1 : middleMiddle - 2*nX + 1;
363 if ((node > -1) && (Map->LID(node) > -1)) {
364 isTouched[i + j*nX + k*nX*nY] = true;
365 numEle += 1;
366 continue;
367 }
368 node = (k == nZ-1) ? -1 : middleMiddle + (2*nX-1)*(2*nY-1) ;
369 if ((node > -1) && (Map->LID(node) > -1)) {
370 isTouched[i + j*nX + k*nX*nY] = true;
371 numEle += 1;
372 continue;
373 }
374 node = (k == 0) ? -1 : middleMiddle - (2*nX-1)*(2*nY-1);
375 if ((node > -1) && (Map->LID(node) > -1)) {
376 isTouched[i + j*nX + k*nX*nY] = true;
377 numEle += 1;
378 continue;
379 }
380 node = ((i == 0) || (j == nY-1)) ? -1 :
381 middleMiddle + 2*nX - 1 - 1;
382 if ((node > -1) && (Map->LID(node) > -1)) {
383 isTouched[i + j*nX + k*nX*nY] = true;
384 numEle += 1;
385 continue;
386 }
387 node = ((i == nX-1) || (j == nY-1)) ? -1 :
388 middleMiddle + 2*nX - 1 + 1;
389 if ((node > -1) && (Map->LID(node) > -1)) {
390 isTouched[i + j*nX + k*nX*nY] = true;
391 numEle += 1;
392 continue;
393 }
394 node = ((i == nX-1) || (j == 0)) ? -1 : middleMiddle - 2*nX + 1 + 1;
395 if ((node > -1) && (Map->LID(node) > -1)) {
396 isTouched[i + j*nX + k*nX*nY] = true;
397 numEle += 1;
398 continue;
399 }
400 node = ((i == 0) || (j == 0)) ? -1 : middleMiddle - 2*nX + 1 - 1;
401 if ((node > -1) && (Map->LID(node) > -1)) {
402 isTouched[i + j*nX + k*nX*nY] = true;
403 numEle += 1;
404 continue;
405 }
406 node = ((i == 0) || (k == nZ-1)) ? -1 : middleMiddle + (2*nX-1)*(2*nY-1) - 1;
407 if ((node > -1) && (Map->LID(node) > -1)) {
408 isTouched[i + j*nX + k*nX*nY] = true;
409 numEle += 1;
410 continue;
411 }
412 node = ((j == nY-1) || (k == nZ-1)) ? -1 : middleMiddle + (2*nX-1)*(2*nY-1) + 2*nX-1;
413 if ((node > -1) && (Map->LID(node) > -1)) {
414 isTouched[i + j*nX + k*nX*nY] = true;
415 numEle += 1;
416 continue;
417 }
418 node = ((i == nX-1) || (k == nZ-1)) ? -1 : middleMiddle + (2*nX-1)*(2*nY-1) + 1;
419 if ((node > -1) && (Map->LID(node) > -1)) {
420 isTouched[i + j*nX + k*nX*nY] = true;
421 numEle += 1;
422 continue;
423 }
424 node = ((j == 0) || (k == nZ-1)) ? -1 : middleMiddle + (2*nX-1)*(2*nY-1) - 2*nX+1;
425 if ((node > -1) && (Map->LID(node) > -1)) {
426 isTouched[i + j*nX + k*nX*nY] = true;
427 numEle += 1;
428 continue;
429 }
430 node = ((i == 0) || (k == 0)) ? -1 : middleMiddle - (2*nX-1)*(2*nY-1) - 1;
431 if ((node > -1) && (Map->LID(node) > -1)) {
432 isTouched[i + j*nX + k*nX*nY] = true;
433 numEle += 1;
434 continue;
435 }
436 node = ((j == nY-1) || (k == 0)) ? -1 : middleMiddle - (2*nX-1)*(2*nY-1) + 2*nX-1;
437 if ((node > -1) && (Map->LID(node) > -1)) {
438 isTouched[i + j*nX + k*nX*nY] = true;
439 numEle += 1;
440 continue;
441 }
442 node = ((i == nX-1) || (k == 0)) ? -1 : middleMiddle - (2*nX-1)*(2*nY-1) + 1;
443 if ((node > -1) && (Map->LID(node) > -1)) {
444 isTouched[i + j*nX + k*nX*nY] = true;
445 numEle += 1;
446 continue;
447 }
448 node = ((j == 0) || (k == 0)) ? -1 : middleMiddle - (2*nX-1)*(2*nY-1) - 2*nX+1;
449 if ((node > -1) && (Map->LID(node) > -1)) {
450 isTouched[i + j*nX + k*nX*nY] = true;
451 numEle += 1;
452 continue;
453 }
454 node = ((i == 0) || (j == nY-1) || (k == nZ-1)) ? -1 :
455 middleMiddle + (2*nX-1)*(2*nY-1) + 2*nX-2;
456 if ((node > -1) && (Map->LID(node) > -1)) {
457 isTouched[i + j*nX + k*nX*nY] = true;
458 numEle += 1;
459 continue;
460 }
461 node = ((i == nX-1) || (j == nY-1) || (k == nZ-1)) ? -1 :
462 middleMiddle + (2*nX-1)*(2*nY-1) + 2*nX;
463 if ((node > -1) && (Map->LID(node) > -1)) {
464 isTouched[i + j*nX + k*nX*nY] = true;
465 numEle += 1;
466 continue;
467 }
468 node = ((i == nX-1) || (j == 0) || (k == nZ-1)) ? -1 :
469 middleMiddle + (2*nX-1)*(2*nY-1) - 2*nX + 2;
470 if ((node > -1) && (Map->LID(node) > -1)) {
471 isTouched[i + j*nX + k*nX*nY] = true;
472 numEle += 1;
473 continue;
474 }
475 node = ((i == 0) || (j == 0) || (k == nZ-1)) ? -1 :
476 middleMiddle + (2*nX-1)*(2*nY-1) - 2*nX;
477 if ((node > -1) && (Map->LID(node) > -1)) {
478 isTouched[i + j*nX + k*nX*nY] = true;
479 numEle += 1;
480 continue;
481 }
482 node = ((i == 0) || (j == nY-1) || (k == 0)) ? -1 :
483 middleMiddle - (2*nX-1)*(2*nY-1) + 2*nX - 2;
484 if ((node > -1) && (Map->LID(node) > -1)) {
485 isTouched[i + j*nX + k*nX*nY] = true;
486 numEle += 1;
487 continue;
488 }
489 node = ((i == nX-1) || (j == nY-1) || (k == 0)) ? -1 :
490 middleMiddle - (2*nX-1)*(2*nY-1) + 2*nX;
491 if ((node > -1) && (Map->LID(node) > -1)) {
492 isTouched[i + j*nX + k*nX*nY] = true;
493 numEle += 1;
494 continue;
495 }
496 node = ((i == nX-1) || (j == 0) || (k == 0)) ? -1 :
497 middleMiddle - (2*nX-1)*(2*nY-1) - 2*nX + 2;
498 if ((node > -1) && (Map->LID(node) > -1)) {
499 isTouched[i + j*nX + k*nX*nY] = true;
500 numEle += 1;
501 continue;
502 }
503 node = ((i == 0) || (j == 0) || (k == 0)) ? -1 :
504 middleMiddle - (2*nX-1)*(2*nY-1) - 2*nX;
505 if ((node > -1) && (Map->LID(node) > -1)) {
506 isTouched[i + j*nX + k*nX*nY] = true;
507 numEle += 1;
508 continue;
509 }
510 }
511 }
512 }
513
514 return numEle;
515
516}
517
518
519void ModeLaplace3DQ2::makeMyElementsTopology(int *elemTopo, bool *isTouched) {
520
521 // Create the element topology for the elements containing nodes for this processor
522 // Note: Put the flag -1 when the node has a Dirichlet boundary condition
523
524 int i, j, k;
525 int numEle = 0;
526
527 for (k=0; k<nZ; ++k) {
528 for (j=0; j<nY; ++j) {
529 for (i=0; i<nX; ++i) {
530 if (isTouched[i + j*nX + k*nX*nY] == false)
531 continue;
532 int middleMiddle = 2*i + 2*j*(2*nX - 1) + 2*k*(2*nX-1)*(2*nY-1);
533 elemTopo[dofEle*numEle+26] = middleMiddle;
534 elemTopo[dofEle*numEle+25] = (i == nX-1) ? -1 : middleMiddle + 1;
535 elemTopo[dofEle*numEle+24] = (i == 0) ? -1 : middleMiddle - 1;
536 elemTopo[dofEle*numEle+23] = (j == nY-1) ? -1 : middleMiddle + 2*nX - 1;
537 elemTopo[dofEle*numEle+22] = (j == 0) ? -1 : middleMiddle - 2*nX + 1;
538 elemTopo[dofEle*numEle+21] = (k == nZ-1) ? -1 : middleMiddle + (2*nX-1)*(2*nY-1) ;
539 elemTopo[dofEle*numEle+20] = (k == 0) ? -1 : middleMiddle - (2*nX-1)*(2*nY-1);
540 elemTopo[dofEle*numEle+19] = ((i == 0) || (j == nY-1)) ? -1 :
541 elemTopo[dofEle*numEle+23] - 1;
542 elemTopo[dofEle*numEle+18] = ((i == nX-1) || (j == nY-1)) ? -1 :
543 elemTopo[dofEle*numEle+23] + 1;
544 elemTopo[dofEle*numEle+17] = ((i == nX-1) || (j == 0)) ? -1 :
545 elemTopo[dofEle*numEle+22] + 1;
546 elemTopo[dofEle*numEle+16] = ((i == 0) || (j == 0)) ? -1 :
547 elemTopo[dofEle*numEle+22] - 1;
548 elemTopo[dofEle*numEle+15] = ((i == 0) || (k == nZ-1)) ? -1 :
549 elemTopo[dofEle*numEle+21] - 1;
550 elemTopo[dofEle*numEle+14] = ((j == nY-1) || (k == nZ-1)) ? -1 :
551 elemTopo[dofEle*numEle+21] + 2*nX - 1;
552 elemTopo[dofEle*numEle+13] = ((i == nX-1) || (k == nZ-1)) ? -1 :
553 elemTopo[dofEle*numEle+21] + 1;
554 elemTopo[dofEle*numEle+12] = ((j == 0) || (k == nZ-1)) ? -1 :
555 elemTopo[dofEle*numEle+21] - 2*nX + 1;
556 elemTopo[dofEle*numEle+11] = ((i == 0) || (k == 0)) ? -1 :
557 elemTopo[dofEle*numEle+20] - 1;
558 elemTopo[dofEle*numEle+10] = ((j == nY-1) || (k == 0)) ? -1 :
559 elemTopo[dofEle*numEle+20] + 2*nX - 1;
560 elemTopo[dofEle*numEle+ 9] = ((i == nX-1) || (k == 0)) ? -1 :
561 elemTopo[dofEle*numEle+20] + 1;
562 elemTopo[dofEle*numEle+ 8] = ((j == 0) || (k == 0)) ? -1 :
563 elemTopo[dofEle*numEle+20] - 2*nX + 1;
564 elemTopo[dofEle*numEle+ 7] = ((i == 0) || (j == nY-1) || (k == nZ-1)) ? -1 :
565 elemTopo[dofEle*numEle+21] + 2*nX - 2;
566 elemTopo[dofEle*numEle+ 6] = ((i == nX-1) || (j == nY-1) || (k == nZ-1)) ? -1 :
567 elemTopo[dofEle*numEle+21] + 2*nX;
568 elemTopo[dofEle*numEle+ 5] = ((i == nX-1) || (j == 0) || (k == nZ-1)) ? -1 :
569 elemTopo[dofEle*numEle+21] - 2*nX + 2;
570 elemTopo[dofEle*numEle+ 4] = ((i == 0) || (j == 0) || (k == nZ-1)) ? -1 :
571 elemTopo[dofEle*numEle+21] - 2*nX;
572 elemTopo[dofEle*numEle+ 3] = ((i == 0) || (j == nY-1) || (k == 0)) ? -1 :
573 elemTopo[dofEle*numEle+20] + 2*nX - 2;
574 elemTopo[dofEle*numEle+ 2] = ((i == nX-1) || (j == nY-1) || (k == 0)) ? -1 :
575 elemTopo[dofEle*numEle+20] + 2*nX;
576 elemTopo[dofEle*numEle+ 1] = ((i == nX-1) || (j == 0) || (k == 0)) ? -1 :
577 elemTopo[dofEle*numEle+20] - 2*nX + 2;
578 elemTopo[dofEle*numEle ] = ((i == 0) || (j == 0) || (k == 0)) ? -1 :
579 elemTopo[dofEle*numEle+20] - 2*nX;
580 numEle += 1;
581 }
582 }
583 }
584
585}
586
587
588void ModeLaplace3DQ2::makeMyConnectivity(int *elemTopo, int numEle, int *connectivity,
589 int *numNz) {
590
591 // This routine creates the connectivity of each managed node
592 // from the element topology.
593
594 int i, j;
595 int localSize = Map->NumMyElements();
596
597 for (i=0; i<localSize; ++i) {
598 numNz[i] = 0;
599 for (j=0; j<maxConnect; ++j) {
600 connectivity[i*maxConnect + j] = -1;
601 }
602 }
603
604 for (i=0; i<numEle; ++i) {
605 for (j=0; j<dofEle; ++j) {
606 if (elemTopo[dofEle*i+j] == -1)
607 continue;
608 int node = Map->LID(elemTopo[dofEle*i+j]);
609 if (node > -1) {
610 int k;
611 for (k=0; k<dofEle; ++k) {
612 int neighbor = elemTopo[dofEle*i+k];
613 if (neighbor == -1)
614 continue;
615 // Check if this neighbor is already stored
616 int l;
617 for (l=0; l<maxConnect; ++l) {
618 if (neighbor == connectivity[node*maxConnect + l])
619 break;
620 if (connectivity[node*maxConnect + l] == -1) {
621 connectivity[node*maxConnect + l] = neighbor;
622 numNz[node] += 1;
623 break;
624 }
625 } // for (l = 0; l < maxConnect; ++l)
626 } // for (k = 0; k < dofEle; ++k)
627 } // if (node > -1)
628 } // for (j = 0; j < dofEle; ++j)
629 } // for (i = 0; i < numEle; ++i)
630
631}
632
633
634void ModeLaplace3DQ2::makeStiffness(int *elemTopo, int numEle, int *connectivity,
635 int *numNz) {
636
637 // Create Epetra_Matrix for stiffness
638 K = new Epetra_CrsMatrix(Copy, *Map, numNz);
639
640 int i;
641 int localSize = Map->NumMyElements();
642
643 double *values = new double[maxConnect];
644 for (i=0; i<maxConnect; ++i)
645 values[i] = 0.0;
646
647 for (i=0; i<localSize; ++i) {
648 int info = K->InsertGlobalValues(Map->GID(i), numNz[i], values, connectivity+maxConnect*i);
649 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
650 "ModeLaplace3DQ2::makeMass(): InsertGlobalValues() returned error code " << info);
651 }
652
653 // Define the elementary matrix
654 double *kel = new double[dofEle*dofEle];
655 makeElementaryStiffness(kel);
656
657 // Assemble the matrix
658 int *indices = new int[dofEle];
659 int numEntries;
660 int j;
661 for (i=0; i<numEle; ++i) {
662 for (j=0; j<dofEle; ++j) {
663 if (elemTopo[dofEle*i + j] == -1)
664 continue;
665 if (Map->LID(elemTopo[dofEle*i+j]) == -1)
666 continue;
667 numEntries = 0;
668 int k;
669 for (k=0; k<dofEle; ++k) {
670 if (elemTopo[dofEle*i+k] == -1)
671 continue;
672 indices[numEntries] = elemTopo[dofEle*i+k];
673 values[numEntries] = kel[dofEle*j + k];
674 numEntries += 1;
675 }
676 int info = K->SumIntoGlobalValues(elemTopo[dofEle*i+j], numEntries, values, indices);
677 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
678 "ModeLaplace3DQ2::makeMass(): SumIntoGlobalValues() returned error code " << info);
679 }
680 }
681
682 delete[] kel;
683 delete[] values;
684 delete[] indices;
685
686 int info = K->FillComplete();
687 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
688 "ModeLaplace3DQ2::makeMass(): SumIntoGlobalValues() returned error code " << info);
689 info = K->OptimizeStorage();
690 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
691 "ModeLaplace3DQ2::makeMass(): SumIntoGlobalValues() returned error code " << info);
692
693}
694
695
696void ModeLaplace3DQ2::makeElementaryStiffness(double *kel) const {
697
698 int i, j, k;
699
700 double hx = Lx/nX;
701 double hy = Ly/nY;
702 double hz = Lz/nZ;
703
704 for (i=0; i<dofEle; ++i) {
705 for (j=0; j<dofEle; ++j) {
706 kel[j + dofEle*i] = 0.0;
707 }
708 }
709
710 double gaussP[3], gaussW[3];
711 gaussP[0] = - sqrt(3.0/5.0); gaussP[1] = 0.0; gaussP[2] = - gaussP[0];
712 gaussW[0] = 5.0/9.0; gaussW[1] = 8.0/9.0; gaussW[2] = 5.0/9.0;
713 double jac = hx*hy*hz/8.0;
714 double *qx = new double[dofEle];
715 double *qy = new double[dofEle];
716 double *qz = new double[dofEle];
717 for (i=0; i<3; ++i) {
718 double xi = gaussP[i];
719 double wi = gaussW[i];
720 for (j=0; j<3; ++j) {
721 double eta = gaussP[j];
722 double wj = gaussW[j];
723 for (k=0; k<3; ++k) {
724 double zeta = gaussP[k];
725 double wk = gaussW[k];
726
727 // Get the shape functions
728
729 qx[ 0] = 2.0/hx*(xi-0.5)*0.5*eta*(eta-1.0)*0.5*zeta*(zeta-1.0);
730 qy[ 0] = 0.5*xi*(xi-1.0)*2.0/hy*(eta-0.5)*0.5*zeta*(zeta-1.0);
731 qz[ 0] = 0.5*xi*(xi-1.0)*0.5*eta*(eta-1.0)*2.0/hz*(zeta-0.5);
732
733 qx[ 1] = 2.0/hx*(xi+0.5)*0.5*eta*(eta-1.0)*0.5*zeta*(zeta-1.0);
734 qy[ 1] = 0.5*xi*(xi+1.0)*2.0/hy*(eta-0.5)*0.5*zeta*(zeta-1.0);
735 qz[ 1] = 0.5*xi*(xi+1.0)*0.5*eta*(eta-1.0)*2.0/hz*(zeta-0.5);
736
737 qx[ 2] = 2.0/hx*(xi+0.5)*0.5*eta*(eta+1.0)*0.5*zeta*(zeta-1.0);
738 qy[ 2] = 0.5*xi*(xi+1.0)*2.0/hy*(eta+0.5)*0.5*zeta*(zeta-1.0);
739 qz[ 2] = 0.5*xi*(xi+1.0)*0.5*eta*(eta+1.0)*2.0/hz*(zeta-0.5);
740
741 qx[ 3] = 2.0/hx*(xi-0.5)*0.5*eta*(eta+1.0)*0.5*zeta*(zeta-1.0);
742 qy[ 3] = 0.5*xi*(xi-1.0)*2.0/hy*(eta+0.5)*0.5*zeta*(zeta-1.0);
743 qz[ 3] = 0.5*xi*(xi-1.0)*0.5*eta*(eta+1.0)*2.0/hz*(zeta-0.5);
744
745 qx[ 4] = 2.0/hx*(xi-0.5)*0.5*eta*(eta-1.0)*0.5*zeta*(zeta+1.0);
746 qy[ 4] = 0.5*xi*(xi-1.0)*2.0/hy*(eta-0.5)*0.5*zeta*(zeta+1.0);
747 qz[ 4] = 0.5*xi*(xi-1.0)*0.5*eta*(eta-1.0)*2.0/hz*(zeta+0.5);
748
749 qx[ 5] = 2.0/hx*(xi+0.5)*0.5*eta*(eta-1.0)*0.5*zeta*(zeta+1.0);
750 qy[ 5] = 0.5*xi*(xi+1.0)*2.0/hy*(eta-0.5)*0.5*zeta*(zeta+1.0);
751 qz[ 5] = 0.5*xi*(xi+1.0)*0.5*eta*(eta-1.0)*2.0/hz*(zeta+0.5);
752
753 qx[ 6] = 2.0/hx*(xi+0.5)*0.5*eta*(eta+1.0)*0.5*zeta*(zeta+1.0);
754 qy[ 6] = 0.5*xi*(xi+1.0)*2.0/hy*(eta+0.5)*0.5*zeta*(zeta+1.0);
755 qz[ 6] = 0.5*xi*(xi+1.0)*0.5*eta*(eta+1.0)*2.0/hz*(zeta+0.5);
756
757 qx[ 7] = 2.0/hx*(xi-0.5)*0.5*eta*(eta+1.0)*0.5*zeta*(zeta+1.0);
758 qy[ 7] = 0.5*xi*(xi-1.0)*2.0/hy*(eta+0.5)*0.5*zeta*(zeta+1.0);
759 qz[ 7] = 0.5*xi*(xi-1.0)*0.5*eta*(eta+1.0)*2.0/hz*(zeta+0.5);
760
761 qx[ 8] = 2.0/hx*(-2.0*xi)*0.5*eta*(eta-1.0)*0.5*zeta*(zeta-1.0);
762 qy[ 8] = (xi+1.0)*(1.0-xi)*2.0/hy*(eta-0.5)*0.5*zeta*(zeta-1.0);
763 qz[ 8] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta-1.0)*2.0/hz*(zeta-0.5);
764
765 qx[ 9] = 2.0/hx*(xi+0.5)*(eta+1.0)*(1.0-eta)*0.5*zeta*(zeta-1.0);
766 qy[ 9] = 0.5*xi*(xi+1.0)*2.0/hy*(-2.0*eta)*0.5*zeta*(zeta-1.0);
767 qz[ 9] = 0.5*xi*(xi+1.0)*(eta+1.0)*(1.0-eta)*2.0/hz*(zeta-0.5);
768
769 qx[10] = 2.0/hx*(-2.0*xi)*0.5*eta*(eta+1.0)*0.5*zeta*(zeta-1.0);
770 qy[10] = (xi+1.0)*(1.0-xi)*2.0/hy*(eta+0.5)*0.5*zeta*(zeta-1.0);
771 qz[10] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta+1.0)*2.0/hz*(zeta-0.5);
772
773 qx[11] = 2.0/hx*(xi-0.5)*(eta+1.0)*(1.0-eta)*0.5*zeta*(zeta-1.0);
774 qy[11] = 0.5*xi*(xi-1.0)*2.0/hy*(-2.0*eta)*0.5*zeta*(zeta-1.0);
775 qz[11] = 0.5*xi*(xi-1.0)*(eta+1.0)*(1.0-eta)*2.0/hz*(zeta-0.5);
776
777 qx[12] = 2.0/hx*(-2.0*xi)*0.5*eta*(eta-1.0)*0.5*zeta*(zeta+1.0);
778 qy[12] = (xi+1.0)*(1.0-xi)*2.0/hy*(eta-0.5)*0.5*zeta*(zeta+1.0);
779 qz[12] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta-1.0)*2.0/hz*(zeta+0.5);
780
781 qx[13] = 2.0/hx*(xi+0.5)*(eta+1.0)*(1.0-eta)*0.5*zeta*(zeta+1.0);
782 qy[13] = 0.5*xi*(xi+1.0)*2.0/hy*(-2.0*eta)*0.5*zeta*(zeta+1.0);
783 qz[13] = 0.5*xi*(xi+1.0)*(eta+1.0)*(1.0-eta)*2.0/hz*(zeta+0.5);
784
785 qx[14] = 2.0/hx*(-2.0*xi)*0.5*eta*(eta+1.0)*0.5*zeta*(zeta+1.0);
786 qy[14] = (xi+1.0)*(1.0-xi)*2.0/hy*(eta+0.5)*0.5*zeta*(zeta+1.0);
787 qz[14] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta+1.0)*2.0/hz*(zeta+0.5);
788
789 qx[15] = 2.0/hx*(xi-0.5)*(eta+1.0)*(1.0-eta)*0.5*zeta*(zeta+1.0);
790 qy[15] = 0.5*xi*(xi-1.0)*2.0/hy*(-2.0*eta)*0.5*zeta*(zeta+1.0);
791 qz[15] = 0.5*xi*(xi-1.0)*(eta+1.0)*(1.0-eta)*2.0/hz*(zeta+0.5);
792
793 qx[16] = 2.0/hx*(xi-0.5)*0.5*eta*(eta-1.0)*(zeta+1.0)*(1.0-zeta);
794 qy[16] = 0.5*xi*(xi-1.0)*2.0/hy*(eta-0.5)*(zeta+1.0)*(1.0-zeta);
795 qz[16] = 0.5*xi*(xi-1.0)*0.5*eta*(eta-1.0)*2.0/hz*(-2.0*zeta);
796
797 qx[17] = 2.0/hx*(xi+0.5)*0.5*eta*(eta-1.0)*(zeta+1.0)*(1.0-zeta);
798 qy[17] = 0.5*xi*(xi+1.0)*2.0/hy*(eta-0.5)*(zeta+1.0)*(1.0-zeta);
799 qz[17] = 0.5*xi*(xi+1.0)*0.5*eta*(eta-1.0)*2.0/hz*(-2.0*zeta);
800
801 qx[18] = 2.0/hx*(xi+0.5)*0.5*eta*(eta+1.0)*(zeta+1.0)*(1.0-zeta);
802 qy[18] = 0.5*xi*(xi+1.0)*2.0/hy*(eta+0.5)*(zeta+1.0)*(1.0-zeta);
803 qz[18] = 0.5*xi*(xi+1.0)*0.5*eta*(eta+1.0)*2.0/hz*(-2.0*zeta);
804
805 qx[19] = 2.0/hx*(xi-0.5)*0.5*eta*(eta+1.0)*(zeta+1.0)*(1.0-zeta);
806 qy[19] = 0.5*xi*(xi-1.0)*2.0/hy*(eta+0.5)*(zeta+1.0)*(1.0-zeta);
807 qz[19] = 0.5*xi*(xi-1.0)*0.5*eta*(eta+1.0)*2.0/hz*(-2.0*zeta);
808
809 qx[20] = 2.0/hx*(-2.0*xi)*(eta+1.0)*(1.0-eta)*0.5*zeta*(zeta-1.0);
810 qy[20] = (xi+1.0)*(1.0-xi)*2.0/hy*(-2.0*eta)*0.5*zeta*(zeta-1.0);
811 qz[20] = (xi+1.0)*(1.0-xi)*(eta+1.0)*(1.0-eta)*2.0/hz*(zeta-0.5);
812
813 qx[21] = 2.0/hx*(-2.0*xi)*(eta+1.0)*(1.0-eta)*0.5*zeta*(zeta+1.0);
814 qy[21] = (xi+1.0)*(1.0-xi)*2.0/hy*(-2.0*eta)*0.5*zeta*(zeta+1.0);
815 qz[21] = (xi+1.0)*(1.0-xi)*(eta+1.0)*(1.0-eta)*2.0/hz*(zeta+0.5);
816
817 qx[22] = 2.0/hx*(-2.0*xi)*0.5*eta*(eta-1.0)*(zeta+1.0)*(1.0-zeta);
818 qy[22] = (xi+1.0)*(1.0-xi)*2.0/hy*(eta-0.5)*(zeta+1.0)*(1.0-zeta);
819 qz[22] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta-1.0)*2.0/hz*(-2.0*zeta);
820
821 qx[23] = 2.0/hx*(-2.0*xi)*0.5*eta*(eta+1.0)*(zeta+1.0)*(1.0-zeta);
822 qy[23] = (xi+1.0)*(1.0-xi)*2.0/hy*(eta+0.5)*(zeta+1.0)*(1.0-zeta);
823 qz[23] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta+1.0)*2.0/hz*(-2.0*zeta);
824
825 qx[24] = 2.0/hx*(xi-0.5)*(eta+1.0)*(1.0-eta)*(zeta+1.0)*(1.0-zeta);
826 qy[24] = 0.5*xi*(xi-1.0)*2.0/hy*(-2.0*eta)*(zeta+1.0)*(1.0-zeta);
827 qz[24] = 0.5*xi*(xi-1.0)*(eta+1.0)*(1.0-eta)*2.0/hz*(-2.0*zeta);
828
829 qx[25] = 2.0/hx*(xi+0.5)*(eta+1.0)*(1.0-eta)*(zeta+1.0)*(1.0-zeta);
830 qy[25] = 0.5*xi*(xi+1.0)*2.0/hy*(-2.0*eta)*(zeta+1.0)*(1.0-zeta);
831 qz[25] = 0.5*xi*(xi+1.0)*(eta+1.0)*(1.0-eta)*2.0/hz*(-2.0*zeta);
832
833 qx[26] = 2.0/hx*(-2.0*xi)*(eta+1.0)*(1.0-eta)*(zeta+1.0)*(1.0-zeta);
834 qy[26] = (xi+1.0)*(1.0-xi)*2.0/hy*(-2.0*eta)*(zeta+1.0)*(1.0-zeta);
835 qz[26] = (xi+1.0)*(1.0-xi)*(eta+1.0)*(1.0-eta)*2.0/hz*(-2.0*zeta);
836
837 // Add in the elementary matrix
838 int ii, jj;
839 for (ii=0; ii<dofEle; ++ii) {
840 for (jj=ii; jj<dofEle; ++jj) {
841 kel[dofEle*ii + jj] += wi*wj*wk*jac*(qx[ii]*qx[jj] + qy[ii]*qy[jj] + qz[ii]*qz[jj]);
842 kel[dofEle*jj + ii] = kel[dofEle*ii + jj];
843 }
844 }
845 }
846 }
847 }
848 delete[] qx;
849 delete[] qy;
850 delete[] qz;
851
852}
853
854
855void ModeLaplace3DQ2::makeMass(int *elemTopo, int numEle, int *connectivity,
856 int *numNz) {
857
858 // Create Epetra_Matrix for mass
859 M = new Epetra_CrsMatrix(Copy, *Map, numNz);
860
861 int i;
862 int localSize = Map->NumMyElements();
863
864 double *values = new double[maxConnect];
865 for (i=0; i<maxConnect; ++i)
866 values[i] = 0.0;
867 for (i=0; i<localSize; ++i) {
868 int info = M->InsertGlobalValues(Map->GID(i), numNz[i], values, connectivity + maxConnect*i);
869 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
870 "ModeLaplace3DQ2::makeMass(): InsertGlobalValues() returned error code " << info);
871 }
872
873 // Define the elementary matrix
874 double *mel = new double[dofEle*dofEle];
875 makeElementaryMass(mel);
876
877 // Assemble the matrix
878 int *indices = new int[dofEle];
879 int numEntries;
880 int j;
881 for (i=0; i<numEle; ++i) {
882 for (j=0; j<dofEle; ++j) {
883 if (elemTopo[dofEle*i + j] == -1)
884 continue;
885 if (Map->LID(elemTopo[dofEle*i+j]) == -1)
886 continue;
887 numEntries = 0;
888 int k;
889 for (k=0; k<dofEle; ++k) {
890 if (elemTopo[dofEle*i+k] == -1)
891 continue;
892 indices[numEntries] = elemTopo[dofEle*i+k];
893 values[numEntries] = mel[dofEle*j + k];
894 numEntries += 1;
895 }
896 int info = M->SumIntoGlobalValues(elemTopo[dofEle*i+j], numEntries, values, indices);
897 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
898 "ModeLaplace3DQ2::makeMass(): SumIntoGlobalValues() returned error code " << info);
899 }
900 }
901
902 delete[] mel;
903 delete[] values;
904 delete[] indices;
905
906 int info = M->FillComplete();
907 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
908 "ModeLaplace3DQ2::makeMass(): SumIntoGlobalValues() returned error code " << info);
909 info = M->OptimizeStorage();
910 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
911 "ModeLaplace3DQ2::makeMass(): SumIntoGlobalValues() returned error code " << info);
912
913}
914
915
916void ModeLaplace3DQ2::makeElementaryMass(double *mel) const {
917
918 int i, j, k;
919
920 double hx = Lx/nX;
921 double hy = Ly/nY;
922 double hz = Lz/nZ;
923
924 for (i=0; i<dofEle; ++i) {
925 for (j=0; j<dofEle; ++j) {
926 mel[j + dofEle*i] = 0.0;
927 }
928 }
929
930 double gaussP[3], gaussW[3];
931 gaussP[0] = - sqrt(3.0/5.0); gaussP[1] = 0.0; gaussP[2] = - gaussP[0];
932 gaussW[0] = 5.0/9.0; gaussW[1] = 8.0/9.0; gaussW[2] = 5.0/9.0;
933 double jac = hx*hy*hz/8.0;
934 double *q = new double[dofEle];
935 for (i=0; i<3; ++i) {
936 double xi = gaussP[i];
937 double wi = gaussW[i];
938 for (j=0; j<3; ++j) {
939 double eta = gaussP[j];
940 double wj = gaussW[j];
941 for (k=0; k<3; ++k) {
942 double zeta = gaussP[k];
943 double wk = gaussW[k];
944 // Get the shape functions
945 q[ 0] = 0.5*xi*(xi-1.0)*0.5*eta*(eta-1.0)*0.5*zeta*(zeta-1.0);
946 q[ 1] = 0.5*xi*(xi+1.0)*0.5*eta*(eta-1.0)*0.5*zeta*(zeta-1.0);
947 q[ 2] = 0.5*xi*(xi+1.0)*0.5*eta*(eta+1.0)*0.5*zeta*(zeta-1.0);
948 q[ 3] = 0.5*xi*(xi-1.0)*0.5*eta*(eta+1.0)*0.5*zeta*(zeta-1.0);
949 q[ 4] = 0.5*xi*(xi-1.0)*0.5*eta*(eta-1.0)*0.5*zeta*(zeta+1.0);
950 q[ 5] = 0.5*xi*(xi+1.0)*0.5*eta*(eta-1.0)*0.5*zeta*(zeta+1.0);
951 q[ 6] = 0.5*xi*(xi+1.0)*0.5*eta*(eta+1.0)*0.5*zeta*(zeta+1.0);
952 q[ 7] = 0.5*xi*(xi-1.0)*0.5*eta*(eta+1.0)*0.5*zeta*(zeta+1.0);
953 q[ 8] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta-1.0)*0.5*zeta*(zeta-1.0);
954 q[ 9] = 0.5*xi*(xi+1.0)*(eta+1.0)*(1.0-eta)*0.5*zeta*(zeta-1.0);
955 q[10] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta+1.0)*0.5*zeta*(zeta-1.0);
956 q[11] = 0.5*xi*(xi-1.0)*(eta+1.0)*(1.0-eta)*0.5*zeta*(zeta-1.0);
957 q[12] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta-1.0)*0.5*zeta*(zeta+1.0);
958 q[13] = 0.5*xi*(xi+1.0)*(eta+1.0)*(1.0-eta)*0.5*zeta*(zeta+1.0);
959 q[14] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta+1.0)*0.5*zeta*(zeta+1.0);
960 q[15] = 0.5*xi*(xi-1.0)*(eta+1.0)*(1.0-eta)*0.5*zeta*(zeta+1.0);
961 q[16] = 0.5*xi*(xi-1.0)*0.5*eta*(eta-1.0)*(zeta+1.0)*(1.0-zeta);
962 q[17] = 0.5*xi*(xi+1.0)*0.5*eta*(eta-1.0)*(zeta+1.0)*(1.0-zeta);
963 q[18] = 0.5*xi*(xi+1.0)*0.5*eta*(eta+1.0)*(zeta+1.0)*(1.0-zeta);
964 q[19] = 0.5*xi*(xi-1.0)*0.5*eta*(eta+1.0)*(zeta+1.0)*(1.0-zeta);
965 q[20] = (xi+1.0)*(1.0-xi)*(eta+1.0)*(1.0-eta)*0.5*zeta*(zeta-1.0);
966 q[21] = (xi+1.0)*(1.0-xi)*(eta+1.0)*(1.0-eta)*0.5*zeta*(zeta+1.0);
967 q[22] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta-1.0)*(zeta+1.0)*(1.0-zeta);
968 q[23] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta+1.0)*(zeta+1.0)*(1.0-zeta);
969 q[24] = 0.5*xi*(xi-1.0)*(eta+1.0)*(1.0-eta)*(zeta+1.0)*(1.0-zeta);
970 q[25] = 0.5*xi*(xi+1.0)*(eta+1.0)*(1.0-eta)*(zeta+1.0)*(1.0-zeta);
971 q[26] = (xi+1.0)*(1.0-xi)*(eta+1.0)*(1.0-eta)*(zeta+1.0)*(1.0-zeta);
972 // Add in the elementary matrix
973 int ii, jj;
974 for (ii=0; ii<dofEle; ++ii) {
975 for (jj=ii; jj<dofEle; ++jj) {
976 mel[dofEle*ii + jj] += wi*wj*wk*jac*q[ii]*q[jj];
977 mel[dofEle*jj + ii] = mel[dofEle*ii + jj];
978 }
979 }
980 }
981 }
982 }
983 delete[] q;
984
985}
986
987
988double ModeLaplace3DQ2::getFirstMassEigenValue() const {
989
990 double hx = Lx/nX;
991 double hy = Ly/nY;
992 double hz = Lz/nZ;
993
994 // Compute the coefficient alphaz
995 double cosk = cos(M_PI*hz/2/Lz);
996 double a = 2.0*cosk;
997 double b = 4.0 + cos(M_PI*hz/Lz);
998 double c = -2.0*cosk;
999 double delta = sqrt(b*b - 4*a*c);
1000 double alphaz = (-b - delta)*0.5/a;
1001
1002 // Compute the coefficient alphay
1003 double cosj = cos(M_PI*hy/2/Ly);
1004 a = 2.0*cosj;
1005 b = 4.0 + cos(M_PI*hy/Ly);
1006 c = -2.0*cosj;
1007 delta = sqrt(b*b - 4*a*c);
1008 double alphay = (-b - delta)*0.5/a;
1009
1010 // Compute the coefficient alphax
1011 double cosi = cos(M_PI*hx/2/Lx);
1012 a = 2.0*cosi;
1013 b = 4.0 + cos(M_PI*hx/Lx);
1014 c = -2.0*cosi;
1015 delta = sqrt(b*b - 4*a*c);
1016 double alphax = (-b - delta)*0.5/a;
1017
1018 double discrete = hx/15.0*(8.0+2*alphax*cosi);
1019 discrete *= hy/15.0*(8.0+2*alphay*cosj);
1020 discrete *= hz/15.0*(8.0+2*alphaz*cosk);
1021
1022 return discrete;
1023
1024}
1025
1026
1027int ModeLaplace3DQ2::eigenCheck(const Epetra_MultiVector &Q, double *lambda,
1028 double *normWeight) const {
1029
1030 int info = 0;
1031 int qc = Q.NumVectors();
1032 int myPid = MyComm.MyPID();
1033
1034 cout.precision(2);
1035 cout.setf(ios::scientific, ios::floatfield);
1036
1037 // Check orthonormality of eigenvectors
1038 double tmp = myVerify.errorOrthonormality(&Q, M);
1039 if (myPid == 0)
1040 cout << " Maximum coefficient in matrix Q^T M Q - I = " << tmp << endl;
1041
1042 // Print out norm of residuals
1043 myVerify.errorEigenResiduals(Q, lambda, K, M, normWeight);
1044
1045 // Check the eigenvalues
1046 int numX = (int) ceil(sqrt(Lx*Lx*lambda[qc-1]/M_PI/M_PI));
1047 numX = (numX > 2*nX) ? 2*nX : numX;
1048 int numY = (int) ceil(sqrt(Ly*Ly*lambda[qc-1]/M_PI/M_PI));
1049 numY = (numY > 2*nY) ? 2*nY : numY;
1050 int numZ = (int) ceil(sqrt(Lz*Lz*lambda[qc-1]/M_PI/M_PI));
1051 numZ = (numZ > 2*nZ) ? 2*nZ : numZ;
1052 int newSize = (numX-1)*(numY-1)*(numZ-1);
1053 double *discrete = new (nothrow) double[2*newSize];
1054 if (discrete == 0) {
1055 return -1;
1056 }
1057 double *continuous = discrete + newSize;
1058
1059 double hx = Lx/nX;
1060 double hy = Ly/nY;
1061 double hz = Lz/nZ;
1062
1063 int i, j, k;
1064 for (k = 1; k < numZ; ++k) {
1065 // Compute the coefficient alphaz
1066 double cosk = cos(k*M_PI*hz/2/Lz);
1067 double a = cosk*(92.0 - 12.0*cos(k*M_PI*hz/Lz));
1068 double b = 48.0 + 32.0*cos(k*M_PI*hz/Lz);
1069 double c = -160.0*cosk;
1070 double delta = sqrt(b*b - 4*a*c);
1071 double alphaz = ((-b - delta)*0.5/a < 0.0) ? (-b + delta)*0.5/a : (-b - delta)*0.5/a;
1072 for (j = 1; j < numY; ++j) {
1073 // Compute the coefficient alphay
1074 double cosj = cos(j*M_PI*hy/2/Ly);
1075 a = cosj*(92.0 - 12.0*cos(j*M_PI*hy/Ly));
1076 b = 48.0 + 32.0*cos(j*M_PI*hy/Ly);
1077 c = -160.0*cosj;
1078 delta = sqrt(b*b - 4*a*c);
1079 double alphay = ((-b - delta)*0.5/a < 0.0) ? (-b + delta)*0.5/a : (-b - delta)*0.5/a;
1080 for (i = 1; i < numX; ++i) {
1081 // Compute the coefficient alphax
1082 double cosi = cos(i*M_PI*hx/2/Lx);
1083 a = cosi*(92.0 - 12.0*cos(i*M_PI*hx/Lx));
1084 b = 48.0 + 32.0*cos(i*M_PI*hx/Lx);
1085 c = -160.0*cosi;
1086 delta = sqrt(b*b - 4*a*c);
1087 double alphax = ((-b - delta)*0.5/a < 0.0) ? (-b + delta)*0.5/a : (-b - delta)*0.5/a;
1088 // Compute the continuous eigenvalue
1089 int pos = i-1 + (j-1)*(numX-1) + (k-1)*(numX-1)*(numY-1);
1090 continuous[pos] = M_PI*M_PI*(i*i/(Lx*Lx) + j*j/(Ly*Ly) + k*k/(Lz*Lz));
1091 // Compute the discrete eigenvalue
1092 discrete[pos] = 240.0*(1.0-alphax*cosi)/((8.0+2*alphax*cosi)*(3.0*hx*hx));
1093 discrete[pos] += 240.0*(1.0-alphay*cosj)/((8.0+2*alphay*cosj)*(3.0*hy*hy));
1094 discrete[pos] += 240.0*(1.0-alphaz*cosk)/((8.0+2*alphaz*cosk)*(3.0*hz*hz));
1095 }
1096 }
1097 }
1098
1099 // Sort the eigenvalues in ascending order
1100 mySort.sortScalars(newSize, continuous);
1101
1102 int *used = new (nothrow) int[newSize];
1103 if (used == 0) {
1104 delete[] discrete;
1105 return -1;
1106 }
1107
1108 mySort.sortScalars(newSize, discrete, used);
1109
1110 int *index = new (nothrow) int[newSize];
1111 if (index == 0) {
1112 delete[] discrete;
1113 delete[] used;
1114 return -1;
1115 }
1116
1117 for (i=0; i<newSize; ++i) {
1118 index[used[i]] = i;
1119 }
1120 delete[] used;
1121
1122 int nMax = myVerify.errorLambda(continuous, discrete, newSize, lambda, qc);
1123
1124 // Define the exact discrete eigenvectors
1125 int localSize = Map->NumMyElements();
1126 double *vQ = new (nothrow) double[(nMax+1)*localSize + nMax];
1127 if (vQ == 0) {
1128 delete[] discrete;
1129 delete[] index;
1130 info = -1;
1131 return info;
1132 }
1133
1134 double *normL2 = vQ + (nMax+1)*localSize;
1135 Epetra_MultiVector Qex(View, *Map, vQ, localSize, nMax);
1136
1137 if ((myPid == 0) && (nMax > 0)) {
1138 cout << endl;
1139 cout << " --- Relative discretization errors for exact eigenvectors ---" << endl;
1140 cout << endl;
1141 cout << " Cont. Values Disc. Values Error H^1 norm L^2 norm\n";
1142 }
1143
1144 for (k=1; k < numZ; ++k) {
1145 for (j=1; j < numY; ++j) {
1146 for (i=1; i < numX; ++i) {
1147 int pos = i-1 + (j-1)*(numX-1) + (k-1)*(numX-1)*(numY-1);
1148 if (index[pos] < nMax) {
1149 int ii;
1150 for (ii=0; ii<localSize; ++ii) {
1151 // Compute the coefficient alphaz
1152 double cosk = cos(k*M_PI*hz/2/Lz);
1153 double a = cosk*(92.0 - 12.0*cos(k*M_PI*hz/Lz));
1154 double b = 48.0 + 32.0*cos(k*M_PI*hz/Lz);
1155 double c = -160.0*cosk;
1156 double delta = sqrt(b*b - 4*a*c);
1157 double alphaz = ((-b - delta)*0.5/a < 0.0) ? (-b + delta)*0.5/a : (-b - delta)*0.5/a;
1158 // Compute the coefficient alphay
1159 double cosj = cos(j*M_PI*hy/2/Ly);
1160 a = cosj*(92.0 - 12.0*cos(j*M_PI*hy/Ly));
1161 b = 48.0 + 32.0*cos(j*M_PI*hy/Ly);
1162 c = -160.0*cosj;
1163 delta = sqrt(b*b - 4*a*c);
1164 double alphay = ((-b - delta)*0.5/a < 0.0) ? (-b + delta)*0.5/a : (-b - delta)*0.5/a;
1165 // Compute the coefficient alphax
1166 double cosi = cos(i*M_PI*hx/2/Lx);
1167 a = cosi*(92.0 - 12.0*cos(i*M_PI*hx/Lx));
1168 b = 48.0 + 32.0*cos(i*M_PI*hx/Lx);
1169 c = -160.0*cosi;
1170 delta = sqrt(b*b - 4*a*c);
1171 double alphax = ((-b - delta)*0.5/a < 0.0) ? (-b + delta)*0.5/a : (-b - delta)*0.5/a;
1172 // Get the value for this eigenvector
1173 double coeff = sin(i*(M_PI/Lx)*x[ii])*sin(j*(M_PI/Ly)*y[ii])*sin(k*(M_PI/Lz)*z[ii]);
1174 if (fabs(x[ii] - floor(x[ii]/hx+0.5)*hx) < 0.25*hx)
1175 coeff *= alphax;
1176 if (fabs(y[ii] - floor(y[ii]/hy+0.5)*hy) < 0.25*hy)
1177 coeff *= alphay;
1178 if (fabs(z[ii] - floor(z[ii]/hz+0.5)*hz) < 0.25*hz)
1179 coeff *= alphaz;
1180 Qex.ReplaceMyValue(ii, index[pos], coeff);
1181 }
1182 // Normalize Qex against the mass matrix
1183 Epetra_MultiVector MQex(View, *Map, vQ + nMax*localSize, localSize, 1);
1184 Epetra_MultiVector Qi(View, Qex, index[pos], 1);
1185 M->Apply(Qi, MQex);
1186 double mnorm = 0.0;
1187 Qi.Dot(MQex, &mnorm);
1188 Qi.Scale(1.0/sqrt(mnorm));
1189 // Compute the L2 norm
1190 Epetra_MultiVector shapeInt(View, *Map, vQ + nMax*localSize, localSize, 1);
1191 for (ii=0; ii<localSize; ++ii) {
1192 double iX, iY, iZ;
1193 if (fabs(x[ii] - floor(x[ii]/hx+0.5)*hx) < 0.25*hx)
1194 iX = 2.0*sin(i*(M_PI/Lx)*x[ii])/(hx*hx*i*(M_PI/Lx)*i*(M_PI/Lx)*i*(M_PI/Lx))*
1195 sqrt(2.0/Lx)*( 3*hx*i*(M_PI/Lx) - 4*sin(i*(M_PI/Lx)*hx) +
1196 cos(i*(M_PI/Lx)*hx)*hx*i*(M_PI/Lx) );
1197 else
1198 iX = 8.0*sin(i*(M_PI/Lx)*x[ii])/(hx*hx*i*(M_PI/Lx)*i*(M_PI/Lx)*i*(M_PI/Lx))*
1199 sqrt(2.0/Lx)*( 2*sin(i*(M_PI/Lx)*0.5*hx) -
1200 cos(i*(M_PI/Lx)*0.5*hx)*hx*i*(M_PI/Lx));
1201 if (fabs(y[ii] - floor(y[ii]/hy+0.5)*hy) < 0.25*hy)
1202 iY = 2.0*sin(j*(M_PI/Ly)*y[ii])/(hy*hy*j*(M_PI/Ly)*j*(M_PI/Ly)*j*(M_PI/Ly))*
1203 sqrt(2.0/Ly)*( 3*hy*j*(M_PI/Ly) - 4*sin(j*(M_PI/Ly)*hy) +
1204 cos(j*(M_PI/Ly)*hy)*hy*j*(M_PI/Ly) );
1205 else
1206 iY = 8.0*sin(j*(M_PI/Ly)*y[ii])/(hy*hy*j*(M_PI/Ly)*j*(M_PI/Ly)*j*(M_PI/Ly))*
1207 sqrt(2.0/Ly)*( 2*sin(j*(M_PI/Ly)*0.5*hy) -
1208 cos(j*(M_PI/Ly)*0.5*hy)*hy*j*(M_PI/Ly));
1209 if (fabs(z[ii] - floor(z[ii]/hz+0.5)*hz) < 0.25*hz)
1210 iZ = 2.0*sin(k*(M_PI/Lz)*z[ii])/(hz*hz*k*(M_PI/Lz)*k*(M_PI/Lz)*k*(M_PI/Lz))*
1211 sqrt(2.0/Lz)*( 3*hz*k*(M_PI/Lz) - 4*sin(k*(M_PI/Lz)*hz) +
1212 cos(k*(M_PI/Lz)*hz)*hz*k*(M_PI/Lz) );
1213 else
1214 iZ = 8.0*sin(k*(M_PI/Lz)*z[ii])/(hz*hz*k*(M_PI/Lz)*k*(M_PI/Lz)*k*(M_PI/Lz))*
1215 sqrt(2.0/Lz)*( 2*sin(k*(M_PI/Lz)*0.5*hz) -
1216 cos(k*(M_PI/Lz)*0.5*hz)*hz*k*(M_PI/Lz));
1217 shapeInt.ReplaceMyValue(ii, 0, iX*iY*iZ);
1218 }
1219 Qi.Dot(shapeInt, normL2 + index[pos]);
1220 } // if index[pos] < nMax)
1221 } // for (i=1; i < numX; ++i)
1222 } // for (j=1; j < numY; ++j)
1223 } // for (k=1; k < numZ; ++k)
1224
1225 if (myPid == 0) {
1226 for (i = 0; i < nMax; ++i) {
1227 double normH1 = continuous[i]*(1.0 - 2.0*normL2[i]) + discrete[i];
1228 normL2[i] = 2.0 - 2.0*normL2[i];
1229 normH1+= normL2[i];
1230 // Print out the result
1231 if (myPid == 0) {
1232 cout << " ";
1233 cout.width(4);
1234 cout << i+1 << ". ";
1235 cout.setf(ios::scientific, ios::floatfield);
1236 cout.precision(8);
1237 cout << continuous[i] << " " << discrete[i] << " ";
1238 cout.precision(3);
1239 cout << fabs(discrete[i] - continuous[i])/continuous[i] << " ";
1240 cout << sqrt(fabs(normH1)/(continuous[i]+1.0)) << " ";
1241 cout << sqrt(fabs(normL2[i])) << endl;
1242 }
1243 } // for (i = 0; i < nMax; ++i)
1244 } // if (myPid == 0)
1245
1246 delete[] discrete;
1247 delete[] index;
1248
1249 // Check the angles between exact discrete eigenvectors and computed
1250
1251 myVerify.errorSubspaces(Q, Qex, M);
1252
1253 delete[] vQ;
1254
1255 return info;
1256
1257}
1258
1259
1260void ModeLaplace3DQ2::memoryInfo() const {
1261
1262 int myPid = MyComm.MyPID();
1263
1264 Epetra_RowMatrix *Mat = dynamic_cast<Epetra_RowMatrix *>(M);
1265 if ((myPid == 0) && (Mat)) {
1266 cout << " Total number of nonzero entries in mass matrix = ";
1267 cout.width(15);
1268 cout << Mat->NumGlobalNonzeros() << endl;
1269 double memSize = Mat->NumGlobalNonzeros()*(sizeof(double) + sizeof(int));
1270 memSize += 2*Mat->NumGlobalRows()*sizeof(int);
1271 cout << " Memory requested for mass matrix per processor = (EST) ";
1272 cout.precision(2);
1273 cout.width(6);
1274 cout.setf(ios::fixed, ios::floatfield);
1275 cout << memSize/1024.0/1024.0/MyComm.NumProc() << " MB " << endl;
1276 cout << endl;
1277 }
1278
1279 Mat = dynamic_cast<Epetra_RowMatrix *>(K);
1280 if ((myPid == 0) && (Mat)) {
1281 cout << " Total number of nonzero entries in stiffness matrix = ";
1282 cout.width(15);
1283 cout << Mat->NumGlobalNonzeros() << endl;
1284 double memSize = Mat->NumGlobalNonzeros()*(sizeof(double) + sizeof(int));
1285 memSize += 2*Mat->NumGlobalRows()*sizeof(int);
1286 cout << " Memory requested for stiffness matrix per processor = (EST) ";
1287 cout.precision(2);
1288 cout.width(6);
1289 cout.setf(ios::fixed, ios::floatfield);
1290 cout << memSize/1024.0/1024.0/MyComm.NumProc() << " MB " << endl;
1291 cout << endl;
1292 }
1293
1294}
1295
1296
1297void ModeLaplace3DQ2::problemInfo() const {
1298
1299 int myPid = MyComm.MyPID();
1300
1301 if (myPid == 0) {
1302 cout.precision(2);
1303 cout.setf(ios::fixed, ios::floatfield);
1304 cout << " --- Problem definition ---\n\n";
1305 cout << " >> Laplace equation in 3D with homogeneous Dirichlet condition\n";
1306 cout << " >> Domain = [0, " << Lx << "] x [0, " << Ly << "] x [0, " << Lz << "]\n";
1307 cout << " >> Orthogonal mesh uniform per direction with Q2 elements (27 nodes)\n";
1308 cout << endl;
1309 cout << " Global size = " << Map->NumGlobalElements() << endl;
1310 cout << endl;
1311 cout << " Number of elements in [0, " << Lx << "] (X-direction): " << nX << endl;
1312 cout << " Number of elements in [0, " << Ly << "] (Y-direction): " << nY << endl;
1313 cout << " Number of elements in [0, " << Lz << "] (Z-direction): " << nZ << endl;
1314 cout << endl;
1315 cout << " Number of interior nodes in the X-direction: " << 2*nX-1 << endl;
1316 cout << " Number of interior nodes in the Y-direction: " << 2*nY-1 << endl;
1317 cout << " Number of interior nodes in the Z-direction: " << 2*nZ-1 << endl;
1318 cout << endl;
1319 }
1320
1321}
1322
1323
Generated on Thu Oct 9 2025 20:35:53 for Anasazi by doxygen 1.9.8