1 /******************************************************************************
2 * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
3 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
4 *
5 * SPDX-License-Identifier: (Apache-2.0 OR MIT)
6 ******************************************************************************/
7
8 /**************************************************************************
9 **************************************************************************
10 * MLI_SFEI Class functions (simplified FEI)
11 **************************************************************************
12 **************************************************************************/
13
14 #include <string.h>
15 #include <stdlib.h>
16 #include <stdio.h>
17 #include "mli_sfei.h"
18
19 /**************************************************************************
20 * constructor
21 *-----------------------------------------------------------------------*/
22
MLI_SFEI(MPI_Comm mpiComm)23 MLI_SFEI::MLI_SFEI(MPI_Comm mpiComm)
24 {
25 mpiComm_ = mpiComm;
26 outputLevel_ = 1;
27 maxElemBlocks_ = 0;
28 nElemBlocks_ = 0;
29 blkNumElems_ = NULL;
30 blkElemNEqns_ = NULL;
31 blkNodeDofs_ = NULL;
32 blkElemEqnLists_ = NULL;
33 blkElemStiffness_ = NULL;
34 // the following variable is added to counter the fact that
35 // the Sandia FEI called addNumElems starting with blkID = 0
36 // while it called loadElemBlock starting with the actual blkID
37 blkIDBase_ = -1;
38 }
39
40 //*************************************************************************
41 // destructor
42 //-------------------------------------------------------------------------
43
~MLI_SFEI()44 MLI_SFEI::~MLI_SFEI()
45 {
46 int i, j;
47 if ( blkElemEqnLists_ != NULL )
48 {
49 for ( i = 0; i < nElemBlocks_; i++ )
50 {
51 for ( j = 0; j < blkNumElems_[i]; j++ )
52 if ( blkElemEqnLists_[i][j] != NULL )
53 delete [] blkElemEqnLists_[i][j];
54 delete [] blkElemEqnLists_[i];
55 }
56 delete [] blkElemEqnLists_;
57 }
58 if ( blkElemStiffness_ != NULL )
59 {
60 for ( i = 0; i < nElemBlocks_; i++ )
61 {
62 for ( j = 0; j < blkNumElems_[i]; j++ )
63 if ( blkElemStiffness_[i][j] != NULL )
64 delete [] blkElemStiffness_[i][j];
65 delete [] blkElemStiffness_[i];
66 }
67 delete [] blkElemStiffness_;
68 }
69 if ( blkNumElems_ != NULL ) delete [] blkNumElems_;
70 if ( blkElemNEqns_ != NULL ) delete [] blkElemNEqns_;
71 if ( blkNodeDofs_ != NULL ) delete [] blkNodeDofs_;
72 }
73
74 //*************************************************************************
75 // set diagnostics output level
76 //-------------------------------------------------------------------------
77
setOutputLevel(int level)78 int MLI_SFEI::setOutputLevel(int level)
79 {
80 if ( level < 0 )
81 {
82 printf("MLI_SFEI::setOutputLevel ERROR - level should be >= 0.\n");
83 return 0;
84 }
85 outputLevel_ = level;
86 return 1;
87 }
88
89 //*************************************************************************
90 // free up stiffness matrices because it is not going to be used any more
91 //-------------------------------------------------------------------------
92
freeStiffnessMatrices()93 int MLI_SFEI::freeStiffnessMatrices()
94 {
95 int i, j;
96 if ( blkElemStiffness_ != NULL )
97 {
98 for ( i = 0; i < nElemBlocks_; i++ )
99 {
100 for ( j = 0; j < blkNumElems_[i]; j++ )
101 if ( blkElemStiffness_[i][j] != NULL )
102 delete [] blkElemStiffness_[i][j];
103 delete [] blkElemStiffness_[i];
104 }
105 delete [] blkElemStiffness_;
106 }
107 blkElemStiffness_ = NULL;
108 blkIDBase_ = -1;
109 return 0;
110 }
111
112 //*************************************************************************
113 // accumulate number of element information
114 //-------------------------------------------------------------------------
115
addNumElems(int elemBlk,int nElems,int nNodesPerElem)116 int MLI_SFEI::addNumElems(int elemBlk, int nElems, int nNodesPerElem)
117 {
118 int iB, *tempBlkNumElems, *tempBlkElemNEqns, *tempBlkNodeDofs;
119
120 if ( elemBlk != nElemBlocks_ && elemBlk != (nElemBlocks_-1) )
121 {
122 printf("MLI_SFEI::addNumElems ERROR : elemBlk %d(%d) invalid\n",
123 elemBlk,nElemBlocks_);
124 return -1;
125 }
126 if ( blkNumElems_ == NULL )
127 {
128 maxElemBlocks_ = 20;
129 nElemBlocks_ = 0;
130 blkNumElems_ = new int[maxElemBlocks_];
131 blkElemNEqns_ = new int[maxElemBlocks_];
132 blkNodeDofs_ = new int[maxElemBlocks_];
133 for ( iB = 0; iB < maxElemBlocks_; iB++ )
134 {
135 blkNumElems_[iB] = 0;
136 blkElemNEqns_[iB] = 0;
137 blkNodeDofs_[iB] = 0;
138 }
139 }
140 if ( elemBlk >= nElemBlocks_ )
141 {
142 if ( nElemBlocks_ >= maxElemBlocks_ )
143 {
144 tempBlkNumElems = blkNumElems_;
145 tempBlkElemNEqns = blkElemNEqns_;
146 tempBlkNodeDofs = blkNodeDofs_;
147 maxElemBlocks_ += 10;
148 blkNumElems_ = new int[maxElemBlocks_];
149 blkElemNEqns_ = new int[maxElemBlocks_];
150 blkNodeDofs_ = new int[maxElemBlocks_];
151 for ( iB = 0; iB < nElemBlocks_; iB++ )
152 {
153 blkNumElems_[iB] = tempBlkNumElems[iB];
154 blkElemNEqns_[iB] = tempBlkElemNEqns[iB];
155 blkNodeDofs_[iB] = tempBlkNodeDofs[iB];
156 }
157 }
158 blkNumElems_[elemBlk] = nElems;
159 blkElemNEqns_[elemBlk] = nNodesPerElem;
160 }
161 else if ( elemBlk >= 0 ) blkNumElems_[elemBlk] += nElems;
162 if ( elemBlk == nElemBlocks_ ) nElemBlocks_++;
163 return 0;
164 }
165
166 //*************************************************************************
167 // initialize the element connectivities
168 //-------------------------------------------------------------------------
169
loadElemBlock(int blkID,int nElems,const int * elemIDs,const double * const * const * stiff,int nEqnsPerElem,const int * const * eqnIndices)170 int MLI_SFEI::loadElemBlock(int blkID, int nElems, const int* elemIDs,
171 const double *const *const *stiff,
172 int nEqnsPerElem, const int *const * eqnIndices)
173 {
174 (void) elemIDs;
175 int iB, iE, iN, iN2, count, currElem, matSize, *nodeList, elemBlk;
176 double *stiffMat;
177
178 if (blkIDBase_ == -1) blkIDBase_ = blkID;
179 elemBlk = blkID - blkIDBase_;
180 if (nElemBlocks_ <= 0) return 0;
181 if (elemBlk < 0 || elemBlk >= nElemBlocks_)
182 {
183 printf("MLI_SFEI::loadElemBlock ERROR : elemBlk %d invalid\n",elemBlk);
184 return -1;
185 }
186 if (blkElemEqnLists_ == NULL)
187 {
188 for (iB = 0; iB < nElemBlocks_; iB++)
189 {
190 if (blkNumElems_[iB] <= 0)
191 {
192 printf("MLI_SFEI::addNumElems ERROR : some elemBlk has 0 elems\n");
193 return -1;
194 }
195 }
196 blkElemEqnLists_ = new int**[nElemBlocks_];
197 blkElemStiffness_ = new double**[nElemBlocks_];
198 for (iB = 0; iB < nElemBlocks_; iB++)
199 {
200 blkElemEqnLists_[iB] = new int*[blkNumElems_[iB]];
201 blkElemStiffness_[iB] = new double*[blkNumElems_[iB]];
202 for (iE = 0; iE < blkNumElems_[iB]; iE++)
203 {
204 blkElemEqnLists_[iB][iE] = NULL;
205 blkElemStiffness_[iB][iE] = NULL;
206 }
207 blkNumElems_[iB] = 0;
208 }
209 }
210 if (nEqnsPerElem != blkElemNEqns_[elemBlk] &&
211 blkElemNEqns_[elemBlk] != 0)
212 blkNodeDofs_[elemBlk] = nEqnsPerElem / blkElemNEqns_[elemBlk];
213
214 blkElemNEqns_[elemBlk] = nEqnsPerElem;
215 currElem = blkNumElems_[elemBlk];
216 matSize = nEqnsPerElem * nEqnsPerElem;
217
218 for (iE = 0; iE < nElems; iE++)
219 {
220 blkElemEqnLists_[elemBlk][currElem] = new int[nEqnsPerElem];
221 nodeList = blkElemEqnLists_[elemBlk][currElem];
222 for (iN = 0; iN < nEqnsPerElem; iN++)
223 nodeList[iN] = eqnIndices[iE][iN];
224 blkElemStiffness_[elemBlk][currElem] = new double[matSize];
225 stiffMat = blkElemStiffness_[elemBlk][currElem];
226 count = 0;
227 for (iN = 0; iN < nEqnsPerElem; iN++)
228 for (iN2 = 0; iN2 < nEqnsPerElem; iN2++)
229 stiffMat[count++] = stiff[iE][iN2][iN];
230 currElem++;
231 }
232 blkNumElems_[elemBlk] = currElem;
233
234 return 0;
235 }
236
237 //*************************************************************************
238 // get block number of elements
239 //-------------------------------------------------------------------------
240
getBlockNumElems(int blkID)241 int MLI_SFEI::getBlockNumElems(int blkID)
242 {
243 if (blkID < 0 || blkID >= nElemBlocks_)
244 {
245 printf("MLI_SFEI::getBlockNumElems ERROR - invalid blkID.\n");
246 return -1;
247 }
248 return blkNumElems_[blkID];
249 }
250
251 //*************************************************************************
252 // get block number of nodes per element
253 //-------------------------------------------------------------------------
254
getBlockElemNEqns(int blkID)255 int MLI_SFEI::getBlockElemNEqns(int blkID)
256 {
257 if (blkID < 0 || blkID >= nElemBlocks_)
258 {
259 printf("MLI_SFEI::getBlockElemNEqns ERROR - invalid blkID.\n");
260 return -1;
261 }
262 return blkElemNEqns_[blkID];
263 }
264
265 //*************************************************************************
266 // get element block nodelists
267 //-------------------------------------------------------------------------
268
getBlockElemEqnLists(int blkID)269 int **MLI_SFEI::getBlockElemEqnLists(int blkID)
270 {
271 if (blkID < 0 || blkID >= nElemBlocks_)
272 {
273 printf("MLI_SFEI::getBlockElemEqnLists ERROR - invalid blkID.\n");
274 return NULL;
275 }
276 return blkElemEqnLists_[blkID];
277 }
278
279 //*************************************************************************
280 // get block element stiffness matrices
281 //-------------------------------------------------------------------------
282
getBlockElemStiffness(int blkID)283 double **MLI_SFEI::getBlockElemStiffness(int blkID)
284 {
285 if (blkID < 0 || blkID >= nElemBlocks_)
286 {
287 printf("MLI_SFEI::getBlockElemStiffness ERROR - invalid blkID.\n");
288 return NULL;
289 }
290 return blkElemStiffness_[blkID];
291 }
292
293