FEI Version of the Day
Loading...
Searching...
No Matches
HexBeamCR.cpp
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#include <fei_macros.hpp>
10
11#include <test_utils/HexBeamCR.hpp>
12
13HexBeamCR::HexBeamCR(int W, int D, int DofPerNode,
14 int decomp, int numProcs, int localProc)
15 : HexBeam(W, D, DofPerNode, decomp, numProcs, localProc)
16{
17 totalNumElems_ = W*W*D;
18 totalNumNodes_ = (W+1)*(W+1)*(D+1) + (W+1)*(W+1)*(2*numProcs_-1);
19
20 numGlobalDOF_ = totalNumNodes_*dofPerNode_;
21
22 numLocalSlices_ = D/numProcs;
23 int remainder = D%numProcs;
24 firstLocalSlice_ = localProc_*numLocalSlices_;
25
26 switch(decomp) {
27 case HexBeamCR::OneD:
28 if (D < numProcs) {
29 fei::console_out() << "HexBeamCR: too many processors." << FEI_ENDL;
30 inErrorState_ = true;
31 break;
32 }
33 if (localProc < remainder) {
34 ++numLocalSlices_;
35 ++firstLocalSlice_;
36 }
37
38 localNumNodes_ = numNodesPerSlice_*(numLocalSlices_+2);
39
40 localNumElems_ = numElemsPerSlice_*numLocalSlices_;
41 numLocalDOF_ = localNumNodes_*dofPerNode_;
42
43 if (localProc > 0) {
44 firstLocalElem_ = localProc*numLocalSlices_*numElemsPerSlice_;
45
46 firstLocalNode_ = localProc*localNumNodes_;
47
48 if (remainder <= localProc && remainder > 0) {
49 firstLocalElem_ += remainder*numElemsPerSlice_;
50 firstLocalNode_ += remainder*numNodesPerSlice_;
51 }
52 }
53
54 break;
55
56 case HexBeamCR::TwoD:
57 case HexBeamCR::ThreeD:
58 default:
59 fei::console_out() << "HexBeamCR: invalid decomp option: " << decomp
60 <<" aborting." << FEI_ENDL;
61 std::abort();
62 }
63
64 localCRslice_ = firstLocalSlice_ + numLocalSlices_/2;
65 numLocalCRs_ = numNodesPerSlice_;
66 if (localProc_ < numProcs_-1) {
67 numLocalCRs_ += numNodesPerSlice_;
68 }
69
70 numNodesPerCR_ = 2;
71}
72
73HexBeamCR::~HexBeamCR()
74{
75}
76
77int HexBeamCR::getElemConnectivity(int elemID, int* nodeIDs)
78{
79 if (elemID < firstLocalElem_ || elemID > firstLocalElem_+localNumElems_) {
80 return(-1);
81 }
82
83 int whichGlobalSlice = elemID/numElemsPerSlice_;
84 int elemX = elemID%W_;
85 int elemY = (elemID%(W_*W_))/W_;
86 //FEI_COUT << "whichGlobalSlice: " << whichGlobalSlice << FEI_ENDL;
87 int firstElemNode = (whichGlobalSlice + localProc_*2)*numNodesPerSlice_
88 + elemY*(W_+1) + elemX;
89
90 if (whichGlobalSlice >= localCRslice_) {
91 firstElemNode += numNodesPerSlice_;
92 }
93
94 nodeIDs[0] = firstElemNode;
95 nodeIDs[1] = firstElemNode+1;
96 nodeIDs[2] = firstElemNode+W_+1;
97 nodeIDs[3] = nodeIDs[2]+1;
98
99 nodeIDs[4] = nodeIDs[0]+numNodesPerSlice_;
100 nodeIDs[5] = nodeIDs[1]+numNodesPerSlice_;
101 nodeIDs[6] = nodeIDs[2]+numNodesPerSlice_;
102 nodeIDs[7] = nodeIDs[3]+numNodesPerSlice_;
103
104 return(0);
105}
106
107int HexBeamCR::getCRNodes(int** nodeIDs)
108{
109 int offset = 0;
110 int firstCRnode = (localCRslice_+localProc_*2)*numNodesPerSlice_;
111 int i;
112 for(i=0; i<numNodesPerSlice_; ++i) {
113 nodeIDs[offset][0] = firstCRnode+i;
114 nodeIDs[offset++][1] = firstCRnode+i+numNodesPerSlice_;
115 }
116
117 if (localProc_ >= numProcs_-1) return(0);
118
119 int nextCRnode = firstLocalNode_ + localNumNodes_ - numNodesPerSlice_;
120 for(i=0; i<numNodesPerSlice_; ++i) {
121 nodeIDs[offset][0] = nextCRnode+i;
122 nodeIDs[offset++][1] = nextCRnode+i+numNodesPerSlice_;
123 }
124
125 return(0);
126}
127
128int HexBeamCR::getElemStiffnessMatrix(int elemID, double* elemMat)
129{
130 if (elemID < firstLocalElem_ || elemID > firstLocalElem_+localNumElems_) {
131 return(-1);
132 }
133
134 int i, len = nodesPerElem_*dofPerNode_*nodesPerElem_*dofPerNode_;
135
136 for(i=0; i<len; ++i) {
137 elemMat[i] = 0.0;
138 }
139
140 //Should set up some semi-realistic stiffness-matrix coefficients here...
141 //For now just use arbitrary numbers and set it up so the matrix won't be
142 //too ill-conditioned. (This is intended for an assembly test more than
143 //a solver test.)
144
145 //Now set the diagonal to 4.0
146 len = nodesPerElem_*dofPerNode_;
147 for(i=0; i<len; ++i) {
148 int offset = i*len+i;
149 elemMat[offset] = 4.0;
150 }
151
152 //Now set some off-diagonals
153 for(i=0; i<len; ++i) {
154 int offset = i*len+i;
155 if (i>1) {
156 elemMat[offset-2] = -0.5;
157 }
158
159 if (i<len-2) {
160 elemMat[offset+2] = -0.5;
161 }
162
163 if (i>3) {
164 elemMat[offset-4] = -0.1;
165 }
166 if (i<len-4) {
167 elemMat[offset+4] = -0.1;
168 }
169 }
170
171 return(0);
172}
173
174int HexBeamCR::getElemLoadVector(int elemID, double* elemVec)
175{
176 if (elemID < firstLocalElem_ || elemID > firstLocalElem_+localNumElems_) {
177 return(-1);
178 }
179
180 int i, len = nodesPerElem_*dofPerNode_;
181 for(i=0; i<len; ++i) {
182 elemVec[i] = 1.0;
183 }
184
185 return(0);
186}
187
188int HexBeamCR::getNumBCNodes()
189{
190 int numBCNodes = (numLocalSlices_+1)*(W_+1);
191 return( numBCNodes );
192}
193
194int HexBeamCR::getBCNodes(int numNodes, int* nodeIDs)
195{
196 if (numNodes != getNumBCNodes()) {
197 return(-1);
198 }
199
200 int firstBCNode = firstLocalNode_ + W_;
201
202 for(int i=0; i<numNodes; ++i) {
203 nodeIDs[i] = firstBCNode + W_+1;
204 }
205
206 return(0);
207}
208
209int HexBeamCR::getBCGammaValues(int numBCDofs, double* gamma)
210{
211 if (numBCDofs != getNumBCNodes()*dofPerNode_) {
212 return(-1);
213 }
214
215 for(int i=0; i<numBCDofs; ++i) {
216 gamma[i] = 2.0;
217 }
218
219 return(0);
220}
221
222int HexBeamCR::getNumSharedNodes()
223{
224 if (numProcs_ < 2) return(0);
225
226 int numSharedNodes = numNodesPerSlice_;
227 if (localProc_ > 0 && localProc_ < numProcs_-1) {
228 numSharedNodes += numNodesPerSlice_;
229 }
230
231 return(numSharedNodes);
232}
233
234int HexBeamCR::getSharedNodes(int numSharedNodes,
235 int*& sharedNodes,
236 int*& numSharingProcsPerNode,
237 int**& sharingProcs)
238{
239 if (numProcs_ < 2) return(0);
240
241 if (numSharedNodes != getNumSharedNodes()) {
242 return(-1);
243 }
244
245 sharedNodes = new int[numSharedNodes];
246 numSharingProcsPerNode = new int[numSharedNodes];
247 sharingProcs = new int*[numSharedNodes];
248 int* sharingProcVals = new int[numSharedNodes];
249 if (sharedNodes == NULL || numSharingProcsPerNode == NULL ||
250 sharingProcs == NULL || sharingProcVals == NULL) {
251 return(-1);
252 }
253
254 int i;
255 for(i=0; i<numSharedNodes; ++i) {
256 numSharingProcsPerNode[i] = 1;
257 sharingProcs[i] = &(sharingProcVals[i]);
258 }
259
260 int firstSharedNode = firstLocalNode_+numNodesPerSlice_*(numLocalSlices_+2);
261 int offset = 0;
262 //FEI_COUT << localProc_ << ": firstSharedNode: " << firstSharedNode << FEI_ENDL;
263 if (localProc_ < numProcs_ - 1) {
264 for(i=0; i<numNodesPerSlice_; ++i) {
265 sharedNodes[offset] = firstSharedNode+i;
266 sharingProcs[offset++][0] = localProc_+1;
267 }
268 }
269
270 firstSharedNode = firstLocalNode_;
271 //FEI_COUT << localProc_ << ":+1 firstSharedNode: " << firstSharedNode << FEI_ENDL;
272 if (localProc_ > 0) {
273 for(i=0; i<numNodesPerSlice_; ++i) {
274 sharedNodes[offset] = firstSharedNode+i;
275 sharingProcs[offset++][0] = localProc_-1;
276 }
277 }
278
279 return(0);
280}
int localProc(MPI_Comm comm)
std::ostream & console_out()
int numProcs(MPI_Comm comm)