Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 
39 const int ModeLaplace3DQ2::dofEle = 27;
40 const int ModeLaplace3DQ2::maxConnect = 125;
41 #ifndef M_PI
42 const double ModeLaplace3DQ2::M_PI = 3.14159265358979323846;
43 #endif
44 
45 
46 ModeLaplace3DQ2::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 
70 ModeLaplace3DQ2::~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 
99 void 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 
162 void 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 
324 int 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 
519 void 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 
588 void 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 
634 void 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 
696 void 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 
855 void 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 
916 void 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 
988 double 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 
1027 int 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 
1260 void 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 
1297 void 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 
Epetra_MultiVector
Epetra_Map
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Epetra_RowMatrix::NumGlobalNonzeros
virtual int NumGlobalNonzeros() const =0
Epetra_Comm
Epetra_CrsMatrix
Epetra_RowMatrix
Epetra_RowMatrix::NumGlobalRows
virtual int NumGlobalRows() const =0
Generated on Fri Jul 25 2025 09:18:08 for Anasazi by   doxygen 1.8.5