FEI Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
PoissonData.cpp
Go to the documentation of this file.
1/*--------------------------------------------------------------------*/
2/* Copyright 2005 Sandia Corporation. */
3/* Under the terms of Contract DE-AC04-94AL85000, there is a */
4/* non-exclusive license for use of this work by or on behalf */
5/* of the U.S. Government. Export of this program may require */
6/* a license from the United States Government. */
7/*--------------------------------------------------------------------*/
8
9//
10//This is a class that will exercise FEI implementations.
11//
12
13#include <fei_macros.hpp>
14#include <fei_defs.h>
15#include <fei_CSVec.hpp>
16
18
20
21#include <cstdlib>
22#include <cmath>
23
24static int int_sqrt(int x) {
25//a not-safe function for taking the sqrt of a square int.
26 return((int)std::ceil(std::sqrt((double)x)));
27}
28
29//==============================================================================
31 int numProcs, int localProc, int outputLevel)
32{
33 //
34 //PoissonData constructor.
35 //
36 //Arguments:
37 //
38 // L: global square size (number-of-elements along side)
39 // numProcs: number of processors participating in this FEI test.
40 // localProc: local processor number.
41 // outputLevel: affects the amount of screen output.
42 //
43
44 L_ = L;
45
46 startElement_ = 0;
48
49 numProcs_ = numProcs;
50 localProc_ = localProc;
51 outputLevel_ = outputLevel;
52
53 check1();
54
55 elem_ = new Poisson_Elem();
56 int err = elem_->allocateInternals(1);
57 err += elem_->allocateLoad(1);
58 err += elem_->allocateStiffness(1);
59 if (err) messageAbort("Allocation error in element.");
60
62 elemIDsAllocated_ = false;
63
64 numFields_ = NULL;
65 fieldIDs_ = NULL;
66
67 elemIDs_ = NULL;
68
70
73 elemSetID_ = 0;
74 elemFormat_ = 0;
75
78
80}
81
82//==============================================================================
84//
85//Destructor -- delete any allocated memory.
86//
87
89
90 if (elemIDsAllocated_) delete [] elemIDs_;
91
93 delete elem_;
94}
95
96//==============================================================================
98//
99//Private function to be called from the constructor, simply makes sure that
100//the constructor's input arguments were reasonable.
101//
102//If they aren't, a message is printed on standard err, and abort() is called.
103//
104 if (L_ <= 0) messageAbort("bar length L <= 0.");
105 if (numProcs_ <= 0) messageAbort("numProcs <= 0.");
106 if (L_%int_sqrt(numProcs_))
107 messageAbort("L must be an integer multiple of sqrt(numProcs).");
108 if (localProc_ < 0) messageAbort("localProc < 0.");
109 if (localProc_ >= numProcs_) messageAbort("localProc >= numProcs.");
110 if (outputLevel_ < 0) messageAbort("outputLevel < 0.");
111}
112
113//==============================================================================
115//
116//Calculate which elements this processor owns. The element domain is a
117//square, and we can assume that sqrt(numProcs_) divides evenly into
118//L_. We're working with a (logically) 2D processor arrangement.
119//Furthermore, the logical processor layout is such that processor 0 is at
120//the bottom left corner of a 2D grid, and a side of the grid is of length
121//sqrt(numProcs_). The element domain is numbered such that element 1 is at
122//the bottom left corner of the square, and element numbers increase from
123//left to right. i.e., element 1 is in position (1,1), element L is in
124//position (1,L), element L+1 is in position (2,1).
125//
126//Use 1-based numbering for the elements and the x- and y- coordinates in
127//the element grid, but use 0-based numbering for processor IDs and the
128//coordinates in the processor grid.
129//
131
133 if (!elemIDs_) messageAbort("ERROR allocating elemIDs_.");
134 elemIDsAllocated_ = true;
135
136 //0-based x-coordinate of this processor in the 2D processor grid.
138
139 //0-based maximum processor x-coordinate.
141
142 //0-based y-coordinate of this processor in the 2D processor grid.
144
145 //0-based maximum processor y-coordinate.
147
148 int sqrtElems = int_sqrt(numLocalElements_);
149 int sqrtProcs = int_sqrt(numProcs_);
150
151 //1-based first-element-on-this-processor
152 startElement_ = 1 + procY_*sqrtProcs*numLocalElements_ + procX_*sqrtElems;
153
154 if (outputLevel_>1) {
155 FEI_COUT << localProc_ << ", calcDist.: numLocalElements: "
156 << numLocalElements_ << ", startElement: " << startElement_
157 << FEI_ENDL;
158 FEI_COUT << localProc_ << ", procX: " << procX_ << ", procY_: " << procY_
159 << ", maxProcX: " << maxProcX_ << ", maxProcY: " << maxProcY_
160 << FEI_ENDL;
161 }
162
163 int offset = 0;
164 for(int i=0; i<sqrtElems; i++) {
165 for(int j=0; j<sqrtElems; j++) {
166 elemIDs_[offset] = (GlobalID)(startElement_ + i*L_ + j);
167 offset++;
168 }
169 }
170}
171
172//==============================================================================
173void PoissonData::messageAbort(const char* message) {
174 fei::console_out() << FEI_ENDL << "PoissonData: " << message
175 << FEI_ENDL << " Aborting." << FEI_ENDL;
176 std::abort();
177}
178
179//==============================================================================
181{
182 //set the elemID on the internal Poisson_Elem instance.
183 elem_->setElemID(elemID);
184 elem_->setElemLength(1.0/L_);
185 elem_->setTotalLength(1.0);
186
187 //now get a pointer to the element's connectivity array and
188 //calculate that connectivity (in place).
189 int size = 0;
190 GlobalID* elemConn = elem_->getElemConnPtr(size);
191 if (size == 0) messageAbort("loadElements: bad conn ptr.");
192
193 calculateConnectivity(elemConn, size, elemID);
194
195 return(elemConn);
196}
197
198//==============================================================================
200{
201 elem_->setElemID(elemID);
202 elem_->setElemLength(1.0/L_);
203 elem_->setTotalLength(1.0);
204
205 //now get a pointer to this element's connectivity array and
206 //calculate that connectivity (in place).
207 int size = 0;
208 GlobalID* elemConn = elem_->getElemConnPtr(size);
209 if (size == 0) messageAbort("loadElemStiffnesses: bad conn ptr.");
210
211 calculateConnectivity(elemConn, size, elemID);
212
214
215 if (outputLevel_>1) {
216 double* x = elem_->getNodalX(size);
217 double* y = elem_->getNodalY(size);
218 FEI_COUT << localProc_ << ", elemID " << elemID << ", nodes: ";
219 for(int j=0; j<size; j++) {
220 FEI_COUT << elemConn[j] << " ";
221 FEI_COUT << "("<<x[j]<<","<<y[j]<<") ";
222 }
224 }
225
227
228 return( elem_->getElemStiff(size) );
229}
230
231//==============================================================================
233{
234 elem_->setElemID(elemID);
235 elem_->setElemLength(1.0/L_);
236 elem_->setTotalLength(1.0);
237
238 //now get a pointer to this element's connectivity array and
239 //calculate that connectivity (in place).
240 int size = 0;
241 GlobalID* elemConn = elem_->getElemConnPtr(size);
242 if (size == 0) messageAbort("loadElemLoads: bad conn ptr.");
243
244 calculateConnectivity(elemConn, size, elemID);
245
247
248 if (outputLevel_>1) {
249 double* x = elem_->getNodalX(size);
250 double* y = elem_->getNodalY(size);
251 FEI_COUT << localProc_ << ", elemID " << elemID << ", nodes: ";
252 for(int j=0; j<size; j++) {
253 FEI_COUT << elemConn[j] << " ";
254 FEI_COUT << "("<<x[j]<<","<<y[j]<<") ";
255 }
257 }
258
260
261 return( elem_->getElemLoad(size));
262
263}
264
265//==============================================================================
267 GlobalID elemID) {
268//
269//Calculate a single element's connectivity array -- the list of nodes
270//that it 'contains'.
271//
272//Note that we're assuming the element is a 2D square.
273//
274 //elemX will be the global 'x-coordinate' of this element in the square. The
275 //'origin' is the lower-left corner of the bar, which is element 1,
276 //and it is in position 1,1.
277 int elemX = (int)elemID%L_;
278 if (elemX == 0) elemX = L_;
279
280 //elemY will be the global (1-based) 'y-coordinate'.
281 int elemY = ((int)elemID - elemX)/L_ + 1;
282
283 //These are the four nodes for this element.
284 GlobalID lowerLeft = elemID + (GlobalID)(elemY-1);
285 GlobalID lowerRight = lowerLeft + (GlobalID)1;
286 GlobalID upperRight = lowerRight + (GlobalID)(L_+1);
287 GlobalID upperLeft = upperRight - (GlobalID)1;
288
289 (void)size;
290
291 //now fill the connectivity array. We'll always fill the connectivity
292 //array with the lower left node first, and then proceed counter-clockwise.
293 conn[0] = lowerLeft;
294 conn[1] = lowerRight;
295 conn[2] = upperRight;
296 conn[3] = upperLeft;
297}
298
299//==============================================================================
301//
302//Set up the field-descriptor variables that will be passed
303//to the FEI's initFields function, beginInitElemBlock function, etc.
304//
305//Note we've hardwired 1 dof per field.
306//
307 fieldSize_ = 1;
308 numFields_ = new int[nodesPerElement_];
309 fieldIDs_ = new int*[nodesPerElement_];
310 for(int i=0; i<nodesPerElement_; i++) {
312 fieldIDs_[i] = new int[fieldsPerNode_];
313 for(int j=0; j<fieldsPerNode_; j++) {
315 }
316 }
318}
319
320//==============================================================================
322
324
325 for(int i=0; i<nodesPerElement_; i++) {
326 delete [] fieldIDs_[i];
327 }
328
329 delete [] fieldIDs_;
330 delete [] numFields_;
331 }
332 fieldArraysAllocated_ = false;
333}
334
335//==============================================================================
336void PoissonData::getLeftSharedNodes(int& numShared, GlobalID* sharedNodeIDs,
337 int* numProcsPerSharedNode,
338 int** sharingProcs) {
339//
340//This function decides whether any of the nodes along the left edge,
341//including the top node but not the bottom node, are shared. It also
342//decides which processors the nodes are shared with.
343//
344
345 if (numProcs_ == 1) {
346 numShared = 0;
347 return;
348 }
349
350 if (procX_ == 0) {
351 //if this proc is on the left edge of the square...
352
353 if (procY_ < maxProcY_) {
354 //if this proc is not the top left proc...
355
356 numShared = 1;
357
358 int topLeftElemIndex = numLocalElements_ -
360
361 elem_->setElemID(elemIDs_[topLeftElemIndex]);
362
363 //now get a pointer to this element's connectivity array and
364 //calculate that connectivity (in place).
365 int size = 0;
366 GlobalID* elemConn = elem_->getElemConnPtr(size);
367 if (size == 0) messageAbort("loadElements: bad conn ptr.");
368
369 calculateConnectivity(elemConn, size, elemIDs_[topLeftElemIndex]);
370
371 sharedNodeIDs[0] = elemConn[3]; //elem's top left node is node 3
372 numProcsPerSharedNode[0] = 2;
373 sharingProcs[0][0] = localProc_;
374 sharingProcs[0][1] = localProc_ + int_sqrt(numProcs_);
375
376 return;
377 }
378 else {
379 //else this proc is the top left proc...
380 numShared = 0;
381 }
382 }
383 else {
384 //else this proc is not on the left edge of the square...
385
386 numShared = int_sqrt(numLocalElements_);
387 int lowerLeftElemIndex = 0;
388
389 int sqrtElems = int_sqrt(numLocalElements_);
390
391 int shOffset = 0;
392 for(int i=0; i<sqrtElems; i++){
393 //stride up the left edge of the local elements...
394 int size=0;
395
396 int elemIndex = lowerLeftElemIndex+i*sqrtElems;
397
398 elem_->setElemID(elemIDs_[elemIndex]);
399
400 //now get a pointer to this element's connectivity array and
401 //calculate that connectivity (in place).
402 GlobalID* nodes = elem_->getElemConnPtr(size);
403 if (size == 0) messageAbort(": bad conn ptr.");
404
405 calculateConnectivity(nodes, size, elemIDs_[elemIndex]);
406
407 //now put in the top left node
408 sharedNodeIDs[shOffset] = nodes[3];
409 sharingProcs[shOffset][0] = localProc_-1;
410 sharingProcs[shOffset][1] = localProc_;
411 numProcsPerSharedNode[shOffset++] = 2;
412 }
413
414 if (procY_ < maxProcY_) {
415 //if this proc isn't on the top edge, the upper left node (the
416 //last one we put into the shared node list) is shared by 4 procs.
417 shOffset--;
418 numProcsPerSharedNode[shOffset] = 4;
419 sharingProcs[shOffset][2] = localProc_ + int_sqrt(numProcs_);
420 sharingProcs[shOffset][3] = sharingProcs[shOffset][2] - 1;
421 }
422 }
423}
424
425//==============================================================================
426void PoissonData::getRightSharedNodes(int& numShared, GlobalID* sharedNodeIDs,
427 int* numProcsPerSharedNode,
428 int** sharingProcs) {
429//
430//This function decides whether any of the nodes along the right edge,
431//including the bottom node but not the top node, are shared. It also
432//decides which processors the nodes are shared with.
433//
434
435 if (numProcs_ == 1) {
436 numShared = 0;
437 return;
438 }
439
440 if (procX_ == maxProcX_) {
441 //if this proc is on the right edge of the square...
442
443 if (procY_ > 0) {
444 //if this proc is not the bottom right proc...
445
446 numShared = 1;
447
448 int lowerRightElemIndex = int_sqrt(numLocalElements_) - 1;
449
450 elem_->setElemID(elemIDs_[lowerRightElemIndex]);
451
452 //now get a pointer to this element's connectivity array and
453 //calculate that connectivity (in place).
454 int size;
455 GlobalID* nodes = elem_->getElemConnPtr(size);
456 if (size == 0) messageAbort(": bad conn ptr.");
457
458 calculateConnectivity(nodes, size, elemIDs_[lowerRightElemIndex]);
459
460 sharedNodeIDs[0] = nodes[1]; //elem's bottom right node is node 1
461 numProcsPerSharedNode[0] = 2;
462 sharingProcs[0][0] = localProc_;
463 sharingProcs[0][1] = localProc_ - int_sqrt(numProcs_);
464
465 return;
466 }
467 else {
468 //else this proc is the bottom right proc...
469 numShared = 0;
470 }
471 }
472 else {
473 //else this proc is not on the right edge of the square...
474
475 numShared = int_sqrt(numLocalElements_);
476 int upperRightElemIndex = numLocalElements_ - 1;
477
478 int sqrtElems = int_sqrt(numLocalElements_);
479
480 int shOffset = 0;
481 for(int i=0; i<sqrtElems; i++){
482 //stride down the right edge of the local elements...
483 int size=0;
484 int elemIndex = upperRightElemIndex-i*sqrtElems;
485 elem_->setElemID(elemIDs_[elemIndex]);
486
487 //now get a pointer to this element's connectivity array and
488 //calculate that connectivity (in place).
489 GlobalID* nodes = elem_->getElemConnPtr(size);
490 if (size == 0) messageAbort(": bad conn ptr.");
491
492 calculateConnectivity(nodes, size, elemIDs_[elemIndex]);
493
494 //now put in the lower right node
495 sharedNodeIDs[shOffset] = nodes[1];
496 sharingProcs[shOffset][0] = localProc_+1;
497 sharingProcs[shOffset][1] = localProc_;
498 numProcsPerSharedNode[shOffset++] = 2;
499 }
500
501 if (procY_ > 0) {
502 //if this proc isn't on the bottom edge, the lower right node (the
503 //last one we put into the shared node list) is shared by 4 procs.
504 shOffset--;
505 numProcsPerSharedNode[shOffset] = 4;
506 sharingProcs[shOffset][2] = localProc_ - int_sqrt(numProcs_);
507 sharingProcs[shOffset][3] = sharingProcs[shOffset][2] + 1;
508 }
509 }
510}
511
512//==============================================================================
513void PoissonData::getTopSharedNodes(int& numShared, GlobalID* sharedNodeIDs,
514 int* numProcsPerSharedNode,
515 int** sharingProcs) {
516//
517//This function decides whether any of the nodes along the top edge,
518//including the right node but not the left node, are shared. It also
519//decides which processors the nodes are shared with.
520//
521
522 if (numProcs_ == 1) {
523 numShared = 0;
524 return;
525 }
526
527 if (procY_ == maxProcY_) {
528 //if this proc is on the top edge of the square...
529
530 if (procX_ < maxProcX_) {
531 //if this proc is not the top right proc...
532
533 numShared = 1;
534
535 int elemIndex = numLocalElements_ - 1;
536
537 elem_->setElemID(elemIDs_[elemIndex]);
538
539 //now get a pointer to this element's connectivity array and
540 //calculate that connectivity (in place).
541 int size;
542 GlobalID* nodes = elem_->getElemConnPtr(size);
543 if (size == 0) messageAbort(": bad conn ptr.");
544
545 calculateConnectivity(nodes, size, elemIDs_[elemIndex]);
546
547 sharedNodeIDs[0] = nodes[2]; //elem's top right node is node 2
548 numProcsPerSharedNode[0] = 2;
549 sharingProcs[0][0] = localProc_;
550 sharingProcs[0][1] = localProc_ + 1;
551
552 return;
553 }
554 else {
555 //else this proc is the top right proc...
556 numShared = 0;
557 }
558 }
559 else {
560 //else this proc is not on the top edge of the square...
561
562 numShared = int_sqrt(numLocalElements_);
563 int topLeftElemIndex = numLocalElements_ - int_sqrt(numLocalElements_);
564
565 int sqrtElems = int_sqrt(numLocalElements_);
566
567 int shOffset = 0;
568 for(int i=0; i<sqrtElems; i++){
569 //stride across the top edge of the local elements...
570 int size=0;
571 int elemIndex = topLeftElemIndex+i;
572
573 elem_->setElemID(elemIDs_[elemIndex]);
574
575 //now get a pointer to this element's connectivity array and
576 //calculate that connectivity (in place).
577 GlobalID* nodes = elem_->getElemConnPtr(size);
578 if (size == 0) messageAbort(": bad conn ptr.");
579
580 calculateConnectivity(nodes, size, elemIDs_[elemIndex]);
581
582 //now put in the upper right node
583 sharedNodeIDs[shOffset] = nodes[2];
584 sharingProcs[shOffset][0] = localProc_+int_sqrt(numProcs_);
585 sharingProcs[shOffset][1] = localProc_;
586 numProcsPerSharedNode[shOffset++] = 2;
587 }
588 if (procX_ < maxProcX_) {
589 //if this proc isn't on the right edge, the top right node (the
590 //last one we put into the shared node list) is shared by 4 procs.
591 shOffset--;
592 numProcsPerSharedNode[shOffset] = 4;
593 sharingProcs[shOffset][2] = localProc_ + 1;
594 sharingProcs[shOffset][3] = sharingProcs[shOffset][0] + 1;
595 }
596 }
597}
598
599//==============================================================================
600void PoissonData::getBottomSharedNodes(int& numShared, GlobalID* sharedNodeIDs,
601 int* numProcsPerSharedNode,
602 int** sharingProcs) {
603//
604//This function decides whether any of the nodes along the bottom edge,
605//including the left node but not the right node, are shared. It also
606//decides which processors the nodes are shared with.
607//
608
609 if (numProcs_ == 1) {
610 numShared = 0;
611 return;
612 }
613
614 if (procY_ == 0) {
615 //if this proc is on the bottom edge of the square...
616
617 if (procX_ > 0) {
618 //if this proc is not the bottom left proc...
619
620 numShared = 1;
621
622 int elemIndex = 0;
623
624 elem_->setElemID(elemIDs_[elemIndex]);
625
626 //now get a pointer to this element's connectivity array and
627 //calculate that connectivity (in place).
628 int size;
629 GlobalID* nodes = elem_->getElemConnPtr(size);
630 if (size == 0) messageAbort(": bad conn ptr.");
631
632 calculateConnectivity(nodes, size, elemIDs_[elemIndex]);
633
634 sharedNodeIDs[0] = nodes[0]; //elem's bottom left node is node 0
635 numProcsPerSharedNode[0] = 2;
636 sharingProcs[0][0] = localProc_;
637 sharingProcs[0][1] = localProc_ - 1;
638
639 return;
640 }
641 else {
642 //else this proc is the top right proc...
643 numShared = 0;
644 }
645 }
646 else {
647 //else this proc is not on the bottom edge of the square...
648
649 numShared = int_sqrt(numLocalElements_);
650 int lowerRightElemIndex = int_sqrt(numLocalElements_) - 1;
651
652 int sqrtElems = int_sqrt(numLocalElements_);
653
654 int shOffset = 0;
655 for(int i=0; i<sqrtElems; i++){
656 //stride across the bottom edge of the local elements, from
657 //right to left...
658 int size=0;
659 int elemIndex = lowerRightElemIndex-i;
660
661 elem_->setElemID(elemIDs_[elemIndex]);
662
663 //now get a pointer to this element's connectivity array and
664 //calculate that connectivity (in place).
665 GlobalID* nodes = elem_->getElemConnPtr(size);
666 if (size == 0) messageAbort(": bad conn ptr.");
667
668 calculateConnectivity(nodes, size, elemIDs_[elemIndex]);
669
670 //now put in the lower left node
671 sharedNodeIDs[shOffset] = nodes[0];
672 sharingProcs[shOffset][0] = localProc_ - int_sqrt(numProcs_);
673 sharingProcs[shOffset][1] = localProc_;
674 numProcsPerSharedNode[shOffset++] = 2;
675 }
676 if (procX_ > 0) {
677 //if this proc isn't on the left edge, the lower left node (the
678 //last one we put into the shared node list) is shared by 4 procs.
679 shOffset--;
680 numProcsPerSharedNode[shOffset] = 4;
681 sharingProcs[shOffset][2] = localProc_ - 1;
682 sharingProcs[shOffset][3] = sharingProcs[shOffset][0] - 1;
683 }
684 }
685}
686
687//==============================================================================
688void PoissonData::printSharedNodes(const char* str,
689 int numShared, GlobalID* nodeIDs,
690 int** shareProcs, int* numShareProcs)
691{
692 for(int i=0; i<numShared; i++) {
693 FEI_COUT << localProc_ << ", " << str << " node: " << (int) nodeIDs[i];
694 FEI_COUT << ", procs: ";
695 for(int j=0; j<numShareProcs[i]; j++) {
696 FEI_COUT << shareProcs[i][j] << " ";
697 }
699 }
700}
701
702//==============================================================================
704//
705//This function figures out which nodes lie on the boundary. The ones that
706//do are added to the BC set, along with appropriate alpha/beta/gamma values.
707//
708 for(int i=0; i<numLocalElements_; i++) {
710 elem_->setElemLength(1.0/L_);
711 elem_->setTotalLength(1.0);
712
713 //now get a pointer to this element's connectivity array and
714 //calculate that connectivity (in place).
715 int size = 0;
716 GlobalID* nodeIDs = elem_->getElemConnPtr(size);
717 if (size == 0) messageAbort("loadElements: bad conn ptr.");
718
719 calculateConnectivity(nodeIDs, size, elemIDs_[i]);
720
722
723 double* xcoord = elem_->getNodalX(size);
724 double* ycoord = elem_->getNodalY(size);
725
726 //now loop over the nodes and see if any are on a boundary.
727 for(int j=0; j<size; j++) {
728 if ((std::abs(xcoord[j]) < 1.e-49) || (std::abs(xcoord[j] - 1.0) < 1.e-49) ||
729 (std::abs(ycoord[j]) < 1.e-49) || (std::abs(ycoord[j] - 1.0) < 1.e-49)) {
730
731 addBCNode(nodeIDs[j], xcoord[j], ycoord[j]);
732 }
733 }
734 }
735}
736
737//==============================================================================
738void PoissonData::addBCNode(GlobalID nodeID, double x, double y){
739
740 std::vector<GlobalID>::iterator
741 iter = std::lower_bound(BCNodeIDs_.begin(), BCNodeIDs_.end(), nodeID);
742
743 if (iter == BCNodeIDs_.end() || *iter != nodeID) {
744 unsigned offset = iter - BCNodeIDs_.begin();
745 BCNodeIDs_.insert(iter, nodeID);
746
747 double bcValue = std::pow(x, 2.0) + std::pow(y, 2.0);
748
749 BCValues_.insert(BCValues_.begin()+offset, bcValue);
750 }
751}
752
753//==============================================================================
755{
756 //first load the information that defines this element block, and
757 //the topology of each element in this element block.
758
759 GlobalID elemBlockID = poissonData.getElemBlockID();
760 int numLocalElements = poissonData.getNumLocalElements();
761 int numNodesPerElement = poissonData.getNumNodesPerElement();
762 int* numFieldsPerNode = poissonData.getNumFieldsPerNodeList();
763 int** fieldIDsTable = poissonData.getNodalFieldIDsTable();
764
765 CHK_ERR( fei->initElemBlock(elemBlockID,
766 numLocalElements,
767 numNodesPerElement,
768 numFieldsPerNode,
769 fieldIDsTable,
770 0, // no element-centered degrees-of-freedom
771 NULL, //null list of elem-dof fieldIDs
773
774 //now let's loop over all of the local elements, giving their
775 //nodal connectivity lists to the FEI.
776
777 GlobalID* elemIDs = poissonData.getLocalElementIDs();
778
779 for(int elem=0; elem<numLocalElements; elem++) {
780 GlobalID* elemConnectivity =
781 poissonData.getElementConnectivity(elemIDs[elem]);
782
783 CHK_ERR( fei->initElem(elemBlockID, elemIDs[elem], elemConnectivity) );
784 }
785
786 return(0);
787}
788
789//==============================================================================
791{
792 int numLocalElements = poissonData.getNumLocalElements();
793 int maxNumSharedNodes = (int)std::sqrt((double)numLocalElements);
794 GlobalID* sharedNodeIDs = new GlobalID[maxNumSharedNodes];
795 int* numProcsPerSharedNode = new int[maxNumSharedNodes];
796 int** sharingProcs = new int*[maxNumSharedNodes];
797 for(int i=0; i<maxNumSharedNodes; i++) sharingProcs[i] = new int[4];
798
799 int numShared;
800
801 //first, get the shared-node data for the left edge of the local block
802
803 poissonData.getLeftSharedNodes(numShared, sharedNodeIDs,
804 numProcsPerSharedNode, sharingProcs);
805
806 CHK_ERR( fei->initSharedNodes(numShared, sharedNodeIDs,
807 numProcsPerSharedNode, sharingProcs));
808
809 //now, get the shared-node data for the right edge of the local block
810
811 poissonData.getRightSharedNodes(numShared, sharedNodeIDs,
812 numProcsPerSharedNode, sharingProcs);
813
814 CHK_ERR( fei->initSharedNodes(numShared, sharedNodeIDs,
815 numProcsPerSharedNode, sharingProcs));
816
817 //now, get the shared-node data for the bottom edge of the local block
818
819 poissonData.getBottomSharedNodes(numShared, sharedNodeIDs,
820 numProcsPerSharedNode, sharingProcs);
821
822 CHK_ERR( fei->initSharedNodes(numShared, sharedNodeIDs,
823 numProcsPerSharedNode, sharingProcs));
824
825 //finally, get the shared-node data for the top edge of the local block
826
827 poissonData.getTopSharedNodes(numShared, sharedNodeIDs,
828 numProcsPerSharedNode, sharingProcs);
829
830 CHK_ERR( fei->initSharedNodes(numShared, sharedNodeIDs,
831 numProcsPerSharedNode, sharingProcs));
832
833 for(int j=0; j<maxNumSharedNodes; j++) delete [] sharingProcs[j];
834 delete [] sharingProcs;
835 delete [] numProcsPerSharedNode;
836 delete [] sharedNodeIDs;
837
838 return(0);
839}
840
841//==============================================================================
843{
844 GlobalID elemBlockID = poissonData.getElemBlockID();
845 int numLocalElements = poissonData.getNumLocalElements();
846 GlobalID* elemIDs = poissonData.getLocalElementIDs();
847
848 for(int elem=0; elem<numLocalElements; elem++) {
849 GlobalID* elemConnectivity =
850 poissonData.getElementConnectivity(elemIDs[elem]);
851 double** elemStiffness = poissonData.getElemStiffness(elemIDs[elem]);
852
853 CHK_ERR( fei->sumInElemMatrix(elemBlockID, elemIDs[elem],
854 elemConnectivity, elemStiffness,
855 poissonData.getElemFormat()));
856
857 double* elemLoad = poissonData.getElemLoad(elemIDs[elem]);
858
859 CHK_ERR( fei->sumInElemRHS(elemBlockID, elemIDs[elem],
860 elemConnectivity, elemLoad));
861 }
862
863 return(0);
864}
865
866//==============================================================================
868{
869 GlobalID elemBlockID = poissonData.getElemBlockID();
870 int numLocalElements = poissonData.getNumLocalElements();
871 GlobalID* elemIDs = poissonData.getLocalElementIDs();
872
873 int numIDs = poissonData.getNumNodesPerElement();
874
875 int* fieldID = poissonData.getFieldIDs();
876
877 fei::CSVec rhs;
878
879 for(int elem=0; elem<numLocalElements; elem++) {
880 GlobalID* elemConnectivity =
881 poissonData.getElementConnectivity(elemIDs[elem]);
882 double** elemStiffness = poissonData.getElemStiffness(elemIDs[elem]);
883
884 CHK_ERR( fei->sumInElemMatrix(elemBlockID, elemIDs[elem],
885 elemConnectivity, elemStiffness,
886 poissonData.getElemFormat()));
887
888 double* elemLoad = poissonData.getElemLoad(elemIDs[elem]);
889
890 for(int i=0; i<numIDs; ++i) {
891 fei::add_entry(rhs, elemConnectivity[i], elemLoad[i]);
892 }
893 }
894
895 fei->loadComplete();
896
897 fei->putIntoRHS(0, *fieldID, rhs.size(),
898 &(rhs.indices()[0]), &(rhs.coefs()[0]));
899
900 return(0);
901}
902
903//==============================================================================
904int load_BC_data(FEI* fei, PoissonData& poissonData)
905{
906 //first, have the data object generate the BC data
907 poissonData.calculateBCs();
908
909 int numBCNodes = poissonData.getNumBCNodes();
910 GlobalID* nodeIDs = poissonData.getBCNodeIDs();
911 int fieldID = poissonData.getBCFieldID();
912 double* values = poissonData.getBCValues();
913
914 std::vector<int> offsets(numBCNodes, 0);
915
916 CHK_ERR( fei->loadNodeBCs(numBCNodes, nodeIDs, fieldID,
917 &offsets[0], values) );
918
919 return(0);
920}
921
922//==============================================================================
924 PoissonData& poissonData)
925{
926 //first load the information that defines this element block, and
927 //the topology of each element in this element block.
928
929 GlobalID elemBlockID = poissonData.getElemBlockID();
930 int numLocalElements = poissonData.getNumLocalElements();
931 int numNodesPerElement = poissonData.getNumNodesPerElement();
932 int** fieldIDsTable = poissonData.getNodalFieldIDsTable();
933
934 int nodeIDType = 0;
935
936 int patternID =
937 matrixGraph->definePattern(numNodesPerElement,
938 nodeIDType, fieldIDsTable[0][0]);
939
940 CHK_ERR( matrixGraph->initConnectivityBlock(elemBlockID,
941 numLocalElements, patternID) );
942
943 //now let's loop over all of the local elements, giving their
944 //nodal connectivity lists to the matrixGraph object.
945
946 GlobalID* elemIDs = poissonData.getLocalElementIDs();
947
948 for(int elem=0; elem<numLocalElements; elem++) {
949 GlobalID* elemConnectivity =
950 poissonData.getElementConnectivity(elemIDs[elem]);
951
952 CHK_ERR( matrixGraph->initConnectivity(elemBlockID, elemIDs[elem],
953 elemConnectivity) );
954 }
955
956 return(0);
957}
958
959//==============================================================================
960int set_shared_nodes(fei::VectorSpace* nodeSpace, PoissonData& poissonData)
961{
962 int numLocalElements = poissonData.getNumLocalElements();
963 int maxNumSharedNodes = (int)std::sqrt((double)numLocalElements);
964 GlobalID* sharedNodeIDs = new GlobalID[maxNumSharedNodes];
965 int* numProcsPerSharedNode = new int[maxNumSharedNodes];
966 int** sharingProcs = new int*[maxNumSharedNodes];
967 for(int i=0; i<maxNumSharedNodes; i++) sharingProcs[i] = new int[4];
968
969 int numShared;
970
971 //first, get the shared-node data for the left edge of the local block
972
973 poissonData.getLeftSharedNodes(numShared, sharedNodeIDs,
974 numProcsPerSharedNode, sharingProcs);
975 int nodeIDType = 0;
976
977 CHK_ERR( nodeSpace->initSharedIDs(numShared, nodeIDType, sharedNodeIDs,
978 numProcsPerSharedNode, sharingProcs));
979
980 //now, get the shared-node data for the right edge of the local block
981
982 poissonData.getRightSharedNodes(numShared, sharedNodeIDs,
983 numProcsPerSharedNode, sharingProcs);
984
985 CHK_ERR( nodeSpace->initSharedIDs(numShared, nodeIDType, sharedNodeIDs,
986 numProcsPerSharedNode, sharingProcs));
987
988 //now, get the shared-node data for the bottom edge of the local block
989
990 poissonData.getBottomSharedNodes(numShared, sharedNodeIDs,
991 numProcsPerSharedNode, sharingProcs);
992
993 CHK_ERR( nodeSpace->initSharedIDs(numShared, nodeIDType, sharedNodeIDs,
994 numProcsPerSharedNode, sharingProcs));
995
996 //finally, get the shared-node data for the top edge of the local block
997
998 poissonData.getTopSharedNodes(numShared, sharedNodeIDs,
999 numProcsPerSharedNode, sharingProcs);
1000
1001 CHK_ERR( nodeSpace->initSharedIDs(numShared, nodeIDType, sharedNodeIDs,
1002 numProcsPerSharedNode, sharingProcs));
1003
1004 for(int j=0; j<maxNumSharedNodes; j++) delete [] sharingProcs[j];
1005 delete [] sharingProcs;
1006 delete [] numProcsPerSharedNode;
1007 delete [] sharedNodeIDs;
1008
1009 return(0);
1010}
1011
1012//==============================================================================
1014 fei::Matrix* mat, fei::Vector* rhs,
1015 PoissonData& poissonData)
1016{
1017 GlobalID elemBlockID = poissonData.getElemBlockID();
1018 int numLocalElements = poissonData.getNumLocalElements();
1019 GlobalID* elemIDs = poissonData.getLocalElementIDs();
1020
1021 int numIndices = matrixGraph->getConnectivityNumIndices(elemBlockID);
1022
1023 std::vector<int> indicesArray(numIndices);
1024 int* indicesPtr = &indicesArray[0];
1025
1026 for(int elem=0; elem<numLocalElements; elem++) {
1027 double** elemStiffness = poissonData.getElemStiffness(elemIDs[elem]);
1028
1029 int checkNumIndices = 0;
1030 CHK_ERR( matrixGraph->getConnectivityIndices(elemBlockID, elemIDs[elem],
1031 numIndices, indicesPtr,
1032 checkNumIndices) );
1033 if (checkNumIndices != numIndices) return(-1);
1034
1035 CHK_ERR( mat->sumIn(elemBlockID, elemIDs[elem],
1036 elemStiffness));
1037
1038 double* elemLoad = poissonData.getElemLoad(elemIDs[elem]);
1039
1040 CHK_ERR( rhs->sumIn(numIndices, indicesPtr, elemLoad));
1041 }
1042
1043 return(0);
1044}
1045
1046//==============================================================================
1048{
1049 //first, have the data object generate the BC data
1050 poissonData.calculateBCs();
1051
1052 int numBCNodes = poissonData.getNumBCNodes();
1053 GlobalID* nodeIDs = poissonData.getBCNodeIDs();
1054 int fieldID = poissonData.getBCFieldID();
1055 double* values = poissonData.getBCValues();
1056
1057 CHK_ERR( linSys->loadEssentialBCs(numBCNodes, nodeIDs, 0, fieldID, 0, values) );
1058
1059 return(0);
1060}
int init_elem_connectivities(FEI *fei, PoissonData &poissonData)
int load_elem_data(FEI *fei, PoissonData &poissonData)
int set_shared_nodes(FEI *fei, PoissonData &poissonData)
int load_BC_data(FEI *fei, PoissonData &poissonData)
int load_elem_data_putrhs(FEI *fei, PoissonData &poissonData)
static int int_sqrt(int x)
Definition: PoissonData.cpp:24
Definition: FEI.hpp:144
int getNumBCNodes()
Definition: PoissonData.hpp:65
void getTopSharedNodes(int &numShared, GlobalID *sharedNodeIDs, int *numProcsPerSharedNode, int **sharingProcs)
GlobalID getElemBlockID()
Definition: PoissonData.hpp:47
void addBCNode(GlobalID nodeID, double x, double y)
bool elemIDsAllocated_
GlobalID * getElementConnectivity(GlobalID elemID)
GlobalID * getBCNodeIDs()
Definition: PoissonData.hpp:66
GlobalID * getLocalElementIDs()
Definition: PoissonData.hpp:50
void calculateConnectivity(GlobalID *conn, int size, GlobalID elemID)
void initializeFieldStuff()
double * getElemLoad(GlobalID elemID)
GlobalID * elemIDs_
int getBCFieldID()
Definition: PoissonData.hpp:67
void check1()
Definition: PoissonData.cpp:97
int nodesPerElement_
void printSharedNodes(const char *str, int numShared, GlobalID *nodeIDs, int **shareProcs, int *numShareProcs)
void getBottomSharedNodes(int &numShared, GlobalID *sharedNodeIDs, int *numProcsPerSharedNode, int **sharingProcs)
int * getFieldIDs()
Definition: PoissonData.hpp:45
void getLeftSharedNodes(int &numShared, GlobalID *sharedNodeIDs, int *numProcsPerSharedNode, int **sharingProcs)
int numLocalElements_
std::vector< double > BCValues_
void messageAbort(const char *message)
int * getNumFieldsPerNodeList()
Definition: PoissonData.hpp:53
void calculateBCs()
void getRightSharedNodes(int &numShared, GlobalID *sharedNodeIDs, int *numProcsPerSharedNode, int **sharingProcs)
int getNumLocalElements()
Definition: PoissonData.hpp:49
int ** fieldIDs_
double ** getElemStiffness(GlobalID elemID)
PoissonData(int L, int numProcs, int localProc, int outputLevel)
Definition: PoissonData.cpp:30
double * getBCValues()
Definition: PoissonData.hpp:68
std::vector< GlobalID > BCNodeIDs_
int getNumNodesPerElement()
Definition: PoissonData.hpp:51
GlobalID elemBlockID_
int getElemFormat()
Definition: PoissonData.hpp:40
int ** getNodalFieldIDsTable()
Definition: PoissonData.hpp:54
void deleteFieldArrays()
void calculateDistribution()
Poisson_Elem * elem_
Definition: PoissonData.hpp:99
bool fieldArraysAllocated_
int * numFields_
void setElemID(GlobalID gNID)
void setElemLength(double len)
double * getNodalX(int &size)
void calculateCoords()
GlobalID * getElemConnPtr(int &size)
void deleteMemory()
void setTotalLength(double len)
double * getNodalY(int &size)
double * getElemLoad(int &size)
double ** getElemStiff(int &size)
void calculateLoad()
int allocateInternals(int DOF)
int allocateLoad(int DOF)
int allocateStiffness(int DOF)
void calculateStiffness()
std::vector< int > & indices()
Definition: fei_CSVec.hpp:31
size_t size() const
Definition: fei_CSVec.hpp:36
std::vector< double > & coefs()
Definition: fei_CSVec.hpp:33
virtual int loadEssentialBCs(int numIDs, const int *IDs, int idType, int fieldID, int offsetIntoField, const double *prescribedValues)
virtual int definePattern(int numIDs, int idType)=0
virtual int initConnectivity(int blockID, int connectivityID, const int *connectedIdentifiers)=0
virtual int getConnectivityNumIndices(int blockID) const =0
virtual int getConnectivityIndices(int blockID, int connectivityID, int indicesAllocLen, int *indices, int &numIndices)=0
virtual int initConnectivityBlock(int blockID, int numConnectivityLists, int patternID, bool diagonal=false)=0
virtual int sumIn(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, int format=0)=0
int initSharedIDs(int numShared, int idType, const int *sharedIDs, const int *numSharingProcsPerID, const int *sharingProcs)
virtual int sumIn(int numValues, const int *indices, const double *values, int vectorIndex=0)=0
#define CHK_ERR(a)
int GlobalID
Definition: fei_defs.h:60
#define FEI_NODE_MAJOR
Definition: fei_defs.h:89
#define FEI_ENDL
#define FEI_COUT
void add_entry(CSVec &vec, int eqn, double coef)
Definition: fei_CSVec.hpp:56
std::ostream & console_out()