Anasazi Version of the Day
Loading...
Searching...
No Matches
ModeLaplace2DQ2.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 "ModeLaplace2DQ2.h"
36#include "Teuchos_Assert.hpp"
37
38
39const int ModeLaplace2DQ2::dofEle = 9;
40const int ModeLaplace2DQ2::maxConnect = 25;
41#ifndef M_PI
42const double ModeLaplace2DQ2::M_PI = 3.14159265358979323846;
43#endif
44
45
46ModeLaplace2DQ2::ModeLaplace2DQ2(const Epetra_Comm &_Comm, double _Lx, int _nX,
47 double _Ly, int _nY)
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 x(0),
59 y(0)
60 {
61
62 preProcess();
63
64}
65
66
67ModeLaplace2DQ2::~ModeLaplace2DQ2() {
68
69 if (Map)
70 delete Map;
71 Map = 0;
72
73 if (K)
74 delete K;
75 K = 0;
76
77 if (M)
78 delete M;
79 M = 0;
80
81 if (x)
82 delete[] x;
83 x = 0;
84
85 if (y)
86 delete[] y;
87 y = 0;
88
89}
90
91
92void ModeLaplace2DQ2::preProcess() {
93
94 // Create the distribution of equations across processors
95 makeMap();
96
97 // Count the number of elements touched by this processor
98 bool *isTouched = new bool[nX*nY];
99 int i, j;
100 for (j = 0; j < nY; ++j)
101 for (i=0; i<nX; ++i)
102 isTouched[i + j*nX] = false;
103
104 int numEle = countElements(isTouched);
105
106 // Create the mesh
107 int *elemTopo = new int[dofEle*numEle];
108 makeMyElementsTopology(elemTopo, isTouched);
109
110 delete[] isTouched;
111
112 // Get the number of nonZeros per row
113 int localSize = Map->NumMyElements();
114 int *numNz = new int[localSize];
115 int *connectivity = new int[localSize*maxConnect];
116 makeMyConnectivity(elemTopo, numEle, connectivity, numNz);
117
118 // Make the stiffness matrix
119 makeStiffness(elemTopo, numEle, connectivity, numNz);
120
121 // Assemble the mass matrix
122 makeMass(elemTopo, numEle, connectivity, numNz);
123
124 // Free some memory
125 delete[] elemTopo;
126 delete[] numNz;
127 delete[] connectivity;
128
129 // Get the geometrical coordinates of the managed nodes
130 double hx = Lx/nX;
131 double hy = Ly/nY;
132
133 x = new double[localSize];
134 y = new double[localSize];
135
136 for (j=0; j<2*nY-1; ++j) {
137 for (i=0; i<2*nX-1; ++i) {
138 int node = i + j*(2*nX-1);
139 if (Map->LID(node) > -1) {
140 x[Map->LID(node)] = (i+1)*hx*0.5;
141 y[Map->LID(node)] = (j+1)*hy*0.5;
142 }
143 }
144 }
145
146}
147
148
149void ModeLaplace2DQ2::makeMap() {
150
151 int numProc = MyComm.NumProc();
152 int globalSize = (2*nX - 1)*(2*nY - 1);
153 assert(globalSize > numProc);
154
155#ifdef _USE_CHACO
156 // Use the partitioner Chaco to distribute the unknowns
157 int *start = new int[globalSize+1];
158 memset(start, 0, (globalSize+1)*sizeof(int));
159
160 int i, j;
161 for (j=0; j<2*nY-1; ++j) {
162 for (i=0; i<2*nX-1; ++i) {
163 int node = i + j*(2*nX-1);
164 int connect = 1;
165 if (i%2 == 1) {
166 connect *= ((i == 1) || (i == 2*nX-3)) ? 4 : 5;
167 }
168 else {
169 connect *= ((i == 0) || (i == 2*nX-2)) ? 2 : 3;
170 }
171 if (j%2 == 1) {
172 connect *= ((j == 1) || (j == 2*nY-3)) ? 4 : 5;
173 }
174 else {
175 connect *= ((j == 0) || (j == 2*nY-2)) ? 2 : 3;
176 }
177 // Don't count the node itself for Chaco
178 start[node+1] = connect - 1;
179 }
180 }
181
182 for (i = 0; i < globalSize; ++i)
183 start[i+1] += start[i];
184
185 int *adjacency = new int[start[globalSize]];
186 memset(adjacency, 0, start[globalSize]*sizeof(int));
187
188 int *elemTopo = new int[dofEle];
189 for (j=0; j<nY; ++j) {
190 for (i=0; i<nX; ++i) {
191 int middleMiddle = 2*i + 2*j*(2*nX - 1);
192 elemTopo[0] = ((i==0) || (j==0)) ? -1 : middleMiddle - 2*nX;
193 elemTopo[1] = ((i==nX-1) || (j==0)) ? -1 : middleMiddle - 2*nX + 2;
194 elemTopo[2] = ((i==nX-1) || (j==nY-1)) ? -1 : middleMiddle + 2*nX;
195 elemTopo[3] = ((i==0) || (j==nY-1)) ? -1 : middleMiddle + 2*nX - 2;
196 elemTopo[4] = (j==0) ? -1 : middleMiddle - 2*nX + 1;
197 elemTopo[5] = (i==nX-1) ? -1 : middleMiddle + 1;
198 elemTopo[6] = (j==nY-1) ? -1 : middleMiddle + 2*nX - 1;
199 elemTopo[7] = (i==0) ? -1 : middleMiddle - 1;
200 elemTopo[8] = middleMiddle;
201 for (int iD = 0; iD < dofEle; ++iD) {
202 if (elemTopo[iD] == -1)
203 continue;
204 for (int jD = 0; jD < dofEle; ++jD) {
205 int neighbor = elemTopo[jD];
206 // Don't count the node itself for Chaco
207 if ((neighbor == -1) || (iD == jD))
208 continue;
209 // Check if this neighbor is already stored
210 for (int l = start[elemTopo[iD]]; l < start[elemTopo[iD]+1]; ++l) {
211 // Note that Chaco uses a Fortran numbering
212 if (adjacency[l] == neighbor + 1)
213 break;
214 if (adjacency[l] == 0) {
215 // Store the node ID
216 // Note that Chaco uses a Fortran numbering
217 adjacency[l] = elemTopo[jD] + 1;
218 break;
219 }
220 } // for (int l = start[elemTopo[iD]]; l < start[elemTopo[iD]+1]; ++l)
221 } // for (int jD = 0; jD < dofEle; ++jD)
222 } // for (int iD = 0; iD < dofEle; ++iD)
223 }
224 }
225 delete[] elemTopo;
226
227 int nDir[3];
228 nDir[0] = numProc;
229 nDir[1] = 0;
230 nDir[2] = 0;
231 short int *partition = new short int[globalSize];
232 // Call Chaco to partition the matrix
233 interface(globalSize, start, adjacency, 0, 0, 0, 0, 0, 0, 0,
234 partition, 1, 1, nDir, 0, 1, 1, 1, 250, 1, 0.001, 7654321L);
235 // Define the Epetra_Map
236 int localSize = 0;
237 int myPid = MyComm.MyPID();
238 for (i = 0; i < globalSize; ++i) {
239 if (partition[i] == myPid)
240 localSize += 1;
241 }
242 int *myRows = new int[localSize];
243 localSize = 0;
244 for (i = 0; i < globalSize; ++i) {
245 if (partition[i] == myPid) {
246 myRows[localSize] = i;
247 localSize +=1 ;
248 }
249 }
250 delete[] partition;
251 delete[] adjacency;
252 delete[] start;
253 Map = new Epetra_Map(globalSize, localSize, myRows, 0, MyComm);
254 delete[] myRows;
255#else
256 // Create a uniform distribution of the unknowns across processors
257 Map = new Epetra_Map(globalSize, 0, MyComm);
258#endif
259
260}
261
262
263int ModeLaplace2DQ2::countElements(bool *isTouched) {
264
265 // This routine counts and flags the elements that contain the nodes
266 // on this processor.
267
268 int i, j;
269 int numEle = 0;
270
271 for (j=0; j<nY; ++j) {
272 for (i=0; i<nX; ++i) {
273 isTouched[i + j*nX] = false;
274 int middleMiddle = 2*i + 2*j*(2*nX - 1);
275 int node;
276 node = ((i==0) || (j==0)) ? -1 : middleMiddle - 2*nX;
277 if ((node > -1) && (Map->LID(node) > -1)) {
278 isTouched[i+j*nX] = true;
279 numEle += 1;
280 continue;
281 }
282 node = ((i==nX-1) || (j==0)) ? -1 : middleMiddle - 2*nX + 2;
283 if ((node > -1) && (Map->LID(node) > -1)) {
284 isTouched[i+j*nX] = true;
285 numEle += 1;
286 continue;
287 }
288 node = ((i==nX-1) || (j==nY-1)) ? -1 : middleMiddle + 2*nX;
289 if ((node > -1) && (Map->LID(node) > -1)) {
290 isTouched[i+j*nX] = true;
291 numEle += 1;
292 continue;
293 }
294 node = ((i==0) || (j==nY-1)) ? -1 : middleMiddle + 2*nX - 2;
295 if ((node > -1) && (Map->LID(node) > -1)) {
296 isTouched[i+j*nX] = true;
297 numEle += 1;
298 continue;
299 }
300 node = (j==0) ? -1 : middleMiddle - 2*nX + 1;
301 if ((node > -1) && (Map->LID(node) > -1)) {
302 isTouched[i+j*nX] = true;
303 numEle += 1;
304 continue;
305 }
306 node = (i==nX-1) ? -1 : middleMiddle + 1;
307 if ((node > -1) && (Map->LID(node) > -1)) {
308 isTouched[i+j*nX] = true;
309 numEle += 1;
310 continue;
311 }
312 node = (j==nY-1) ? -1 : middleMiddle + 2*nX - 1;
313 if ((node > -1) && (Map->LID(node) > -1)) {
314 isTouched[i+j*nX] = true;
315 numEle += 1;
316 continue;
317 }
318 node = (i==0) ? -1 : middleMiddle - 1;
319 if ((node > -1) && (Map->LID(node) > -1)) {
320 isTouched[i+j*nX] = true;
321 numEle += 1;
322 continue;
323 }
324 node = middleMiddle;
325 if ((node > -1) && (Map->LID(node) > -1)) {
326 isTouched[i+j*nX] = true;
327 numEle += 1;
328 continue;
329 }
330 }
331 }
332
333 return numEle;
334
335}
336
337
338void ModeLaplace2DQ2::makeMyElementsTopology(int *elemTopo, bool *isTouched) {
339
340 // Create the element topology for the elements containing nodes for this processor
341 // Note: Put the flag -1 when the node has a Dirichlet boundary condition
342
343 int i, j;
344 int numEle = 0;
345
346 for (j=0; j<nY; ++j) {
347 for (i=0; i<nX; ++i) {
348 if (isTouched[i + j*nX] == false)
349 continue;
350 int middleMiddle = 2*i + 2*j*(2*nX - 1);
351 elemTopo[dofEle*numEle] = ((i==0) || (j==0)) ? -1 : middleMiddle - 2*nX;
352 elemTopo[dofEle*numEle+1] = ((i==nX-1) || (j==0)) ? -1 : middleMiddle - 2*nX + 2;
353 elemTopo[dofEle*numEle+2] = ((i==nX-1) || (j==nY-1)) ? -1 : middleMiddle + 2*nX;
354 elemTopo[dofEle*numEle+3] = ((i==0) || (j==nY-1)) ? -1 : middleMiddle + 2*nX - 2;
355 elemTopo[dofEle*numEle+4] = (j==0) ? -1 : middleMiddle - 2*nX + 1;
356 elemTopo[dofEle*numEle+5] = (i==nX-1) ? -1 : middleMiddle + 1;
357 elemTopo[dofEle*numEle+6] = (j==nY-1) ? -1 : middleMiddle + 2*nX - 1;
358 elemTopo[dofEle*numEle+7] = (i==0) ? -1 : middleMiddle - 1;
359 elemTopo[dofEle*numEle+8] = middleMiddle;
360 numEle += 1;
361 }
362 }
363
364}
365
366
367void ModeLaplace2DQ2::makeMyConnectivity(int *elemTopo, int numEle, int *connectivity,
368 int *numNz) {
369
370 // This routine creates the connectivity of each managed node
371 // from the element topology.
372
373 int i, j;
374 int localSize = Map->NumMyElements();
375
376 for (i=0; i<localSize; ++i) {
377 numNz[i] = 0;
378 for (j=0; j<maxConnect; ++j) {
379 connectivity[i*maxConnect + j] = -1;
380 }
381 }
382
383 for (i=0; i<numEle; ++i) {
384 for (j=0; j<dofEle; ++j) {
385 if (elemTopo[dofEle*i+j] == -1)
386 continue;
387 int node = Map->LID(elemTopo[dofEle*i+j]);
388 if (node > -1) {
389 int k;
390 for (k=0; k<dofEle; ++k) {
391 int neighbor = elemTopo[dofEle*i+k];
392 if (neighbor == -1)
393 continue;
394 // Check if this neighbor is already stored
395 int l;
396 for (l=0; l<maxConnect; ++l) {
397 if (neighbor == connectivity[node*maxConnect + l])
398 break;
399 if (connectivity[node*maxConnect + l] == -1) {
400 connectivity[node*maxConnect + l] = neighbor;
401 numNz[node] += 1;
402 break;
403 }
404 } // for (l = 0; l < maxConnect; ++l)
405 } // for (k = 0; k < dofEle; ++k)
406 } // if (node > -1)
407 } // for (j = 0; j < dofEle; ++j)
408 } // for (i = 0; i < numEle; ++i)
409
410}
411
412
413void ModeLaplace2DQ2::makeStiffness(int *elemTopo, int numEle, int *connectivity,
414 int *numNz) {
415
416 // Create Epetra_Matrix for stiffness
417 K = new Epetra_CrsMatrix(Copy, *Map, numNz);
418
419 int i;
420 int localSize = Map->NumMyElements();
421
422 double *values = new double[maxConnect];
423 for (i=0; i<maxConnect; ++i)
424 values[i] = 0.0;
425
426 for (i=0; i<localSize; ++i) {
427 int info = K->InsertGlobalValues(Map->GID(i), numNz[i], values, connectivity+maxConnect*i);
428 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
429 "ModeLaplace2DQ2::makeStiffness(): InsertGlobalValues() returned error code " << info);
430 }
431
432 // Define the elementary matrix
433 double *kel = new double[dofEle*dofEle];
434 makeElementaryStiffness(kel);
435
436 // Assemble the matrix
437 int *indices = new int[dofEle];
438 int numEntries;
439 int j;
440 for (i=0; i<numEle; ++i) {
441 for (j=0; j<dofEle; ++j) {
442 if (elemTopo[dofEle*i + j] == -1)
443 continue;
444 if (Map->LID(elemTopo[dofEle*i+j]) == -1)
445 continue;
446 numEntries = 0;
447 int k;
448 for (k=0; k<dofEle; ++k) {
449 if (elemTopo[dofEle*i+k] == -1)
450 continue;
451 indices[numEntries] = elemTopo[dofEle*i+k];
452 values[numEntries] = kel[dofEle*j + k];
453 numEntries += 1;
454 }
455 int info = K->SumIntoGlobalValues(elemTopo[dofEle*i+j], numEntries, values, indices);
456 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
457 "ModeLaplace2DQ2::makeStiffness(): SumIntoGlobalValues() returned error code " << info);
458 }
459 }
460
461 delete[] kel;
462 delete[] values;
463 delete[] indices;
464
465 int info = K->FillComplete();
466 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
467 "ModeLaplace2DQ2::makeStiffness(): FillComplete() returned error code " << info);
468 info = K->OptimizeStorage();
469 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
470 "ModeLaplace2DQ2::makeStiffness(): OptimizeStorage() returned error code " << info);
471
472}
473
474
475void ModeLaplace2DQ2::makeElementaryStiffness(double *kel) const {
476
477 int i, j;
478
479 double hx = Lx/nX;
480 double hy = Ly/nY;
481
482 for (i=0; i<dofEle; ++i) {
483 for (j=0; j<dofEle; ++j) {
484 kel[j + dofEle*i] = 0.0;
485 }
486 }
487
488 double gaussP[3], gaussW[3];
489 gaussP[0] = - sqrt(3.0/5.0); gaussP[1] = 0.0; gaussP[2] = - gaussP[0];
490 gaussW[0] = 5.0/9.0; gaussW[1] = 8.0/9.0; gaussW[2] = 5.0/9.0;
491 double jac = hx*hy/4.0;
492 double *qx = new double[dofEle];
493 double *qy = new double[dofEle];
494 for (i=0; i<3; ++i) {
495 double xi = gaussP[i];
496 double wi = gaussW[i];
497 for (j=0; j<3; ++j) {
498 double eta = gaussP[j];
499 double wj = gaussW[j];
500 // Get the shape functions
501 qx[0] = 2.0/hx*(xi-0.5)*0.5*eta*(eta-1.0);
502 qx[1] = 2.0/hx*(xi+0.5)*0.5*eta*(eta-1.0);
503 qx[2] = 2.0/hx*(xi+0.5)*0.5*eta*(eta+1.0);
504 qx[3] = 2.0/hx*(xi-0.5)*0.5*eta*(eta+1.0);
505 qx[4] = 2.0/hx*(-2.0*xi)*0.5*eta*(eta-1.0);
506 qx[5] = 2.0/hx*(xi+0.5)*(eta+1.0)*(1.0-eta);
507 qx[6] = 2.0/hx*(-2.0*xi)*0.5*eta*(eta+1.0);
508 qx[7] = 2.0/hx*(xi-0.5)*(eta+1.0)*(1.0-eta);
509 qx[8] = 2.0/hx*(-2.0*xi)*(eta+1.0)*(1.0-eta);
510 qy[0] = 2.0/hy*0.5*xi*(xi-1.0)*(eta-0.5);
511 qy[1] = 2.0/hy*0.5*xi*(xi+1.0)*(eta-0.5);
512 qy[2] = 2.0/hy*0.5*xi*(xi+1.0)*(eta+0.5);
513 qy[3] = 2.0/hy*0.5*xi*(xi-1.0)*(eta+0.5);
514 qy[4] = 2.0/hy*(xi+1.0)*(1.0-xi)*(eta-0.5);
515 qy[5] = 2.0/hy*0.5*xi*(xi+1.0)*(-2.0*eta);
516 qy[6] = 2.0/hy*(xi+1.0)*(1.0-xi)*(eta+0.5);
517 qy[7] = 2.0/hy*0.5*xi*(xi-1.0)*(-2.0*eta);
518 qy[8] = 2.0/hy*(xi+1.0)*(1.0-xi)*(-2.0*eta);
519 // Add in the elementary matrix
520 int ii, jj;
521 for (ii=0; ii<dofEle; ++ii) {
522 for (jj=ii; jj<dofEle; ++jj) {
523 kel[dofEle*ii + jj] += wi*wj*jac*(qx[ii]*qx[jj] + qy[ii]*qy[jj]);
524 kel[dofEle*jj + ii] = kel[dofEle*ii + jj];
525 }
526 }
527 }
528 }
529 delete[] qx;
530 delete[] qy;
531
532}
533
534
535void ModeLaplace2DQ2::makeMass(int *elemTopo, int numEle, int *connectivity,
536 int *numNz) {
537
538 // Create Epetra_Matrix for mass
539 M = new Epetra_CrsMatrix(Copy, *Map, numNz);
540
541 int i;
542 int localSize = Map->NumMyElements();
543
544 double *values = new double[maxConnect];
545 for (i=0; i<maxConnect; ++i)
546 values[i] = 0.0;
547 for (i=0; i<localSize; ++i) {
548 int info = M->InsertGlobalValues(Map->GID(i), numNz[i], values, connectivity + maxConnect*i);
549 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
550 "ModeLaplace2DQ2::makeMass(): InsertGlobalValues() returned error code " << info);
551 }
552
553 // Define the elementary matrix
554 double *mel = new double[dofEle*dofEle];
555 makeElementaryMass(mel);
556
557 // Assemble the matrix
558 int *indices = new int[dofEle];
559 int numEntries;
560 int j;
561 for (i=0; i<numEle; ++i) {
562 for (j=0; j<dofEle; ++j) {
563 if (elemTopo[dofEle*i + j] == -1)
564 continue;
565 if (Map->LID(elemTopo[dofEle*i+j]) == -1)
566 continue;
567 numEntries = 0;
568 int k;
569 for (k=0; k<dofEle; ++k) {
570 if (elemTopo[dofEle*i+k] == -1)
571 continue;
572 indices[numEntries] = elemTopo[dofEle*i+k];
573 values[numEntries] = mel[dofEle*j + k];
574 numEntries += 1;
575 }
576 int info = M->SumIntoGlobalValues(elemTopo[dofEle*i+j], numEntries, values, indices);
577 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
578 "ModeLaplace2DQ2::makeMass(): SumIntoGlobalValues() returned error code " << info);
579 }
580 }
581
582 delete[] mel;
583 delete[] values;
584 delete[] indices;
585
586 int info = M->FillComplete();
587 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
588 "ModeLaplace2DQ2::makeMass(): FillComplete() returned error code " << info);
589 info = M->OptimizeStorage();
590 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
591 "ModeLaplace2DQ2::makeMass(): OptimizeStorage() returned error code " << info);
592
593}
594
595
596void ModeLaplace2DQ2::makeElementaryMass(double *mel) const {
597
598 int i, j;
599
600 double hx = Lx/nX;
601 double hy = Ly/nY;
602
603 for (i=0; i<dofEle; ++i) {
604 for (j=0; j<dofEle; ++j) {
605 mel[j + dofEle*i] = 0.0;
606 }
607 }
608
609 double gaussP[3], gaussW[3];
610 gaussP[0] = - sqrt(3.0/5.0); gaussP[1] = 0.0; gaussP[2] = - gaussP[0];
611 gaussW[0] = 5.0/9.0; gaussW[1] = 8.0/9.0; gaussW[2] = 5.0/9.0;
612 double jac = hx*hy/4.0;
613 double *q = new double[dofEle];
614 for (i=0; i<3; ++i) {
615 double xi = gaussP[i];
616 double wi = gaussW[i];
617 for (j=0; j<3; ++j) {
618 double eta = gaussP[j];
619 double wj = gaussW[j];
620 // Get the shape functions
621 q[0] = 0.5*xi*(xi-1.0)*0.5*eta*(eta-1.0);
622 q[1] = 0.5*xi*(xi+1.0)*0.5*eta*(eta-1.0);
623 q[2] = 0.5*xi*(xi+1.0)*0.5*eta*(eta+1.0);
624 q[3] = 0.5*xi*(xi-1.0)*0.5*eta*(eta+1.0);
625 q[4] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta-1.0);
626 q[5] = 0.5*xi*(xi+1.0)*(eta+1.0)*(1.0-eta);
627 q[6] = (xi+1.0)*(1.0-xi)*0.5*eta*(eta+1.0);
628 q[7] = 0.5*xi*(xi-1.0)*(eta+1.0)*(1.0-eta);
629 q[8] = (xi+1.0)*(1.0-xi)*(eta+1.0)*(1.0-eta);
630 // Add in the elementary matrix
631 int ii, jj;
632 for (ii=0; ii<dofEle; ++ii) {
633 for (jj=ii; jj<dofEle; ++jj) {
634 mel[dofEle*ii + jj] += wi*wj*jac*q[ii]*q[jj];
635 mel[dofEle*jj + ii] = mel[dofEle*ii + jj];
636 }
637 }
638 }
639 }
640 delete[] q;
641
642}
643
644
645double ModeLaplace2DQ2::getFirstMassEigenValue() const {
646
647 double hx = Lx/nX;
648 double hy = Ly/nY;
649
650 // Compute the coefficient alphay
651 double cosj = cos(M_PI*hy/2/Ly);
652 double a = 2.0*cosj;
653 double b = 4.0 + cos(M_PI*hy/Ly);
654 double c = -2.0*cosj;
655 double delta = sqrt(b*b - 4*a*c);
656 double alphay = (-b - delta)*0.5/a;
657
658 // Compute the coefficient alphax
659 double cosi = cos(M_PI*hx/2/Lx);
660 a = 2.0*cosi;
661 b = 4.0 + cos(M_PI*hx/Lx);
662 c = -2.0*cosi;
663 delta = sqrt(b*b - 4*a*c);
664 double alphax = (-b - delta)*0.5/a;
665
666 double discrete = hx/15.0*(8.0+2*alphax*cosi);
667 discrete *= hy/15.0*(8.0+2*alphay*cosj);
668
669 return discrete;
670
671}
672
673
674int ModeLaplace2DQ2::eigenCheck(const Epetra_MultiVector &Q, double *lambda,
675 double *normWeight) const {
676
677 int info = 0;
678 int qc = Q.NumVectors();
679 int myPid = MyComm.MyPID();
680
681 cout.precision(2);
682 cout.setf(ios::scientific, ios::floatfield);
683
684 // Check orthonormality of eigenvectors
685 double tmp = myVerify.errorOrthonormality(&Q, M);
686 if (myPid == 0)
687 cout << " Maximum coefficient in matrix Q^T M Q - I = " << tmp << endl;
688
689 // Print out norm of residuals
690 myVerify.errorEigenResiduals(Q, lambda, K, M, normWeight);
691
692 // Check the eigenvalues
693 int numX = (int) ceil(sqrt(Lx*Lx*lambda[qc-1]/M_PI/M_PI));
694 numX = (numX > 2*nX) ? 2*nX : numX;
695 int numY = (int) ceil(sqrt(Ly*Ly*lambda[qc-1]/M_PI/M_PI));
696 numY = (numY > 2*nY) ? 2*nY : numY;
697 int newSize = (numX-1)*(numY-1);
698 double *discrete = new (nothrow) double[2*newSize];
699 if (discrete == 0) {
700 return -1;
701 }
702 double *continuous = discrete + newSize;
703
704 double hx = Lx/nX;
705 double hy = Ly/nY;
706
707 int i, j;
708 for (j = 1; j< numY; ++j) {
709 // Compute the coefficient alphay
710 double cosj = cos(j*M_PI*hy/2/Ly);
711 double a = cosj*(92.0 - 12.0*cos(j*M_PI*hy/Ly));
712 double b = 48.0 + 32.0*cos(j*M_PI*hy/Ly);
713 double c = -160.0*cosj;
714 double delta = sqrt(b*b - 4*a*c);
715 double alphay = ((-b - delta)*0.5/a < 0.0) ? (-b + delta)*0.5/a : (-b - delta)*0.5/a;
716 for (i = 1; i < numX; ++i) {
717 int pos = i-1 + (j-1)*(numX-1);
718 continuous[pos] = M_PI*M_PI*(i*i/(Lx*Lx) + j*j/(Ly*Ly));
719 // Compute the coefficient alphax
720 double cosi = cos(i*M_PI*hx/2/Lx);
721 a = cosi*(92.0 - 12.0*cos(i*M_PI*hx/Lx));
722 b = 48.0 + 32.0*cos(i*M_PI*hx/Lx);
723 c = -160.0*cosi;
724 delta = sqrt(b*b - 4*a*c);
725 double alphax = ((-b - delta)*0.5/a < 0.0) ? (-b + delta)*0.5/a : (-b - delta)*0.5/a;
726 // Compute the discrete eigenvalue
727 discrete[pos] = 240.0*(1.0-alphax*cosi)/((8.0+2*alphax*cosi)*(3.0*hx*hx));
728 discrete[pos] += 240.0*(1.0-alphay*cosj)/((8.0+2*alphay*cosj)*(3.0*hy*hy));
729 }
730 }
731
732 // Sort the eigenvalues in ascending order
733 mySort.sortScalars(newSize, continuous);
734
735 int *used = new (nothrow) int[newSize];
736 if (used == 0) {
737 delete[] discrete;
738 return -1;
739 }
740
741 mySort.sortScalars(newSize, discrete, used);
742
743 int *index = new (nothrow) int[newSize];
744 if (index == 0) {
745 delete[] discrete;
746 delete[] used;
747 return -1;
748 }
749
750 for (i=0; i<newSize; ++i) {
751 index[used[i]] = i;
752 }
753 delete[] used;
754
755 int nMax = myVerify.errorLambda(continuous, discrete, newSize, lambda, qc);
756
757 // Define the exact discrete eigenvectors
758 int localSize = Map->NumMyElements();
759 double *vQ = new (nothrow) double[(nMax+1)*localSize + nMax];
760 if (vQ == 0) {
761 delete[] discrete;
762 delete[] index;
763 info = -1;
764 return info;
765 }
766
767 double *normL2 = vQ + (nMax+1)*localSize;
768 Epetra_MultiVector Qex(View, *Map, vQ, localSize, nMax);
769
770 if ((myPid == 0) && (nMax > 0)) {
771 cout << endl;
772 cout << " --- Relative discretization errors for exact eigenvectors ---" << endl;
773 cout << endl;
774 cout << " Cont. Values Disc. Values Error H^1 norm L^2 norm\n";
775 }
776
777 for (j=1; j < numY; ++j) {
778 for (i=1; i < numX; ++i) {
779 if (index[i-1 + (j-1)*(numX-1)] < nMax) {
780 // Compute the coefficients alphax and alphay
781 double cosj = cos(j*M_PI*hy/2/Ly);
782 double a = cosj*(92.0 - 12.0*cos(j*M_PI*hy/Ly));
783 double b = 48.0 + 32.0*cos(j*M_PI*hy/Ly);
784 double c = -160.0*cosj;
785 double delta = sqrt(b*b - 4*a*c);
786 double alphay = ((-b - delta)*0.5/a < 0.0) ? (-b + delta)*0.5/a : (-b - delta)*0.5/a;
787 double cosi = cos(i*M_PI*hx/2/Lx);
788 a = cosi*(92.0 - 12.0*cos(i*M_PI*hx/Lx));
789 b = 48.0 + 32.0*cos(i*M_PI*hx/Lx);
790 c = -160.0*cosi;
791 delta = sqrt(b*b - 4*a*c);
792 double alphax = ((-b - delta)*0.5/a < 0.0) ? (-b + delta)*0.5/a : (-b - delta)*0.5/a;
793 int ii;
794 for (ii=0; ii<localSize; ++ii) {
795 double coeff = sin(i*(M_PI/Lx)*x[ii])*sin(j*(M_PI/Ly)*y[ii]);
796 if (fabs(x[ii] - floor(x[ii]/hx+0.5)*hx) < 0.25*hx)
797 coeff *= alphax;
798 if (fabs(y[ii] - floor(y[ii]/hy+0.5)*hy) < 0.25*hy)
799 coeff *= alphay;
800 Qex.ReplaceMyValue(ii, index[i-1+(j-1)*(numX-1)], coeff);
801 }
802 // Normalize Qex against the mass matrix
803 Epetra_MultiVector MQex(View, *Map, vQ + nMax*localSize, localSize, 1);
804 Epetra_MultiVector Qi(View, Qex, index[i-1+(j-1)*(numX-1)], 1);
805 M->Apply(Qi, MQex);
806 double mnorm = 0.0;
807 Qi.Dot(MQex, &mnorm);
808 Qi.Scale(1.0/sqrt(mnorm));
809 // Compute the L2 norm
810 Epetra_MultiVector shapeInt(View, *Map, vQ + nMax*localSize, localSize, 1);
811 for (ii=0; ii<localSize; ++ii) {
812 double iX, iY;
813 if (fabs(x[ii] - floor(x[ii]/hx+0.5)*hx) < 0.25*hx)
814 iX = 2.0*sin(i*(M_PI/Lx)*x[ii])/(hx*hx*i*(M_PI/Lx)*i*(M_PI/Lx)*i*(M_PI/Lx))*
815 sqrt(2.0/Lx)*( 3*hx*i*(M_PI/Lx) - 4*sin(i*(M_PI/Lx)*hx) +
816 cos(i*(M_PI/Lx)*hx)*hx*i*(M_PI/Lx) );
817 else
818 iX = 8.0*sin(i*(M_PI/Lx)*x[ii])/(hx*hx*i*(M_PI/Lx)*i*(M_PI/Lx)*i*(M_PI/Lx))*
819 sqrt(2.0/Lx)*( 2*sin(i*(M_PI/Lx)*0.5*hx) -
820 cos(i*(M_PI/Lx)*0.5*hx)*hx*i*(M_PI/Lx));
821 if (fabs(y[ii] - floor(y[ii]/hy+0.5)*hy) < 0.25*hy)
822 iY = 2.0*sin(j*(M_PI/Ly)*y[ii])/(hy*hy*j*(M_PI/Ly)*j*(M_PI/Ly)*j*(M_PI/Ly))*
823 sqrt(2.0/Ly)*( 3*hy*j*(M_PI/Ly) - 4*sin(j*(M_PI/Ly)*hy) +
824 cos(j*(M_PI/Ly)*hy)*hy*j*(M_PI/Ly) );
825 else
826 iY = 8.0*sin(j*(M_PI/Ly)*y[ii])/(hy*hy*j*(M_PI/Ly)*j*(M_PI/Ly)*j*(M_PI/Ly))*
827 sqrt(2.0/Ly)*( 2*sin(j*(M_PI/Ly)*0.5*hy) -
828 cos(j*(M_PI/Ly)*0.5*hy)*hy*j*(M_PI/Ly));
829 shapeInt.ReplaceMyValue(ii, 0, iX*iY);
830 }
831 Qi.Dot(shapeInt, normL2 + index[i-1+(j-1)*(numX-1)]);
832 } // if index[i-1+(j-1)*(numX-1)] < nMax)
833 } // for (i=1; i < numX; ++i)
834 } // for (j=1; j < numY; ++j)
835
836 if (myPid == 0) {
837 for (i = 0; i < nMax; ++i) {
838 double normH1 = continuous[i]*(1.0 - 2.0*normL2[i]) + discrete[i];
839 normL2[i] = 2.0 - 2.0*normL2[i];
840 normH1+= normL2[i];
841 // Print out the result
842 if (myPid == 0) {
843 cout << " ";
844 cout.width(4);
845 cout << i+1 << ". ";
846 cout.setf(ios::scientific, ios::floatfield);
847 cout.precision(8);
848 cout << continuous[i] << " " << discrete[i] << " ";
849 cout.precision(3);
850 cout << fabs(discrete[i] - continuous[i])/continuous[i] << " ";
851 cout << sqrt(fabs(normH1)/(continuous[i]+1.0)) << " ";
852 cout << sqrt(fabs(normL2[i])) << endl;
853 }
854 } // for (i = 0; i < nMax; ++i)
855 } // if (myPid == 0)
856
857 delete[] discrete;
858 delete[] index;
859
860 // Check the angles between exact discrete eigenvectors and computed
861
862 myVerify.errorSubspaces(Q, Qex, M);
863
864 delete[] vQ;
865
866 return info;
867
868}
869
870
871void ModeLaplace2DQ2::memoryInfo() const {
872
873 int myPid = MyComm.MyPID();
874
875 Epetra_RowMatrix *Mat = dynamic_cast<Epetra_RowMatrix *>(M);
876 if ((myPid == 0) && (Mat)) {
877 cout << " Total number of nonzero entries in mass matrix = ";
878 cout.width(15);
879 cout << Mat->NumGlobalNonzeros() << endl;
880 double memSize = Mat->NumGlobalNonzeros()*(sizeof(double) + sizeof(int));
881 memSize += 2*Mat->NumGlobalRows()*sizeof(int);
882 cout << " Memory requested for mass matrix per processor = (EST) ";
883 cout.precision(2);
884 cout.width(6);
885 cout.setf(ios::fixed, ios::floatfield);
886 cout << memSize/1024.0/1024.0/MyComm.NumProc() << " MB " << endl;
887 cout << endl;
888 }
889
890 Mat = dynamic_cast<Epetra_RowMatrix *>(K);
891 if ((myPid == 0) && (Mat)) {
892 cout << " Total number of nonzero entries in stiffness matrix = ";
893 cout.width(15);
894 cout << Mat->NumGlobalNonzeros() << endl;
895 double memSize = Mat->NumGlobalNonzeros()*(sizeof(double) + sizeof(int));
896 memSize += 2*Mat->NumGlobalRows()*sizeof(int);
897 cout << " Memory requested for stiffness matrix per processor = (EST) ";
898 cout.precision(2);
899 cout.width(6);
900 cout.setf(ios::fixed, ios::floatfield);
901 cout << memSize/1024.0/1024.0/MyComm.NumProc() << " MB " << endl;
902 cout << endl;
903 }
904
905}
906
907
908void ModeLaplace2DQ2::problemInfo() const {
909
910 int myPid = MyComm.MyPID();
911
912 if (myPid == 0) {
913 cout.precision(2);
914 cout.setf(ios::fixed, ios::floatfield);
915 cout << " --- Problem definition ---\n\n";
916 cout << " >> Laplace equation in 2D with homogeneous Dirichlet condition\n";
917 cout << " >> Domain = [0, " << Lx << "] x [0, " << Ly << "]\n";
918 cout << " >> Orthogonal mesh uniform per direction with Q2 elements (9 nodes)\n";
919 cout << endl;
920 cout << " Global size = " << Map->NumGlobalElements() << endl;
921 cout << endl;
922 cout << " Number of elements in [0, " << Lx << "] (X-direction): " << nX << endl;
923 cout << " Number of elements in [0, " << Ly << "] (Y-direction): " << nY << endl;
924 cout << endl;
925 cout << " Number of interior nodes in the X-direction: " << 2*nX-1 << endl;
926 cout << " Number of interior nodes in the Y-direction: " << 2*nY-1 << endl;
927 cout << endl;
928 }
929
930}
931
932
Generated on Thu Oct 9 2025 20:35:52 for Anasazi by doxygen 1.9.8