FEI Version of the Day
Loading...
Searching...
No Matches
fei_Filter.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_sstream.hpp"
10
11#include "fei_CommUtils.hpp"
12#include "fei_TemplateUtils.hpp"
13
14#include "fei_defs.h"
15#include "fei_NodeDescriptor.hpp"
16#include "fei_NodeDatabase.hpp"
17#include "fei_BlockDescriptor.hpp"
18#include "SNL_FEI_Structure.hpp"
19#include "snl_fei_Utils.hpp"
20#include "fei_Filter.hpp"
21
22#include <cmath>
23#include <algorithm>
24
25#undef fei_file
26#define fei_file "fei_Filter.cpp"
27
28#include "fei_ErrMacros.hpp"
29
30//------------------------------------------------------------------------------
32 : problemStructure_(probStruct),
33 logInput_(false),
34 logInputStream_(NULL),
35 outputLevel_(0)
36{
37}
38
39//------------------------------------------------------------------------------
42
43//------------------------------------------------------------------------------
44void Filter::setLogStream(FEI_OSTREAM* logstrm)
45{
46 logInputStream_ = logstrm;
47}
48
49//------------------------------------------------------------------------------
50FEI_OSTREAM* Filter::logStream()
51{
52 return( logInputStream_ );
53}
54
55//------------------------------------------------------------------------------
56void Filter::copyStiffness(const double* const* elemStiff,
57 int numRows, int elemFormat,
58 double** copy)
59{
60 //
61 //Unpack the element stiffness array in elemStiff into a full dense
62 //'copy'.
63 //
64 int i, j;
65 const double* elStiff_i = NULL;
66
67 switch (elemFormat) {
68
69 case FEI_DENSE_ROW:
70 for (i = 0; i < numRows; i++) {
71 elStiff_i = elemStiff[i];
72 for (j = 0; j < numRows; j++) {
73 copy[i][j] = elStiff_i[j];
74 }
75 }
76 break;
77
78 case FEI_UPPER_SYMM_ROW:
79 for (i = 0; i < numRows; i++) {
80 elStiff_i = elemStiff[i];
81 int jcol=0;
82 for (j = i; j < numRows; j++) {
83 copy[i][j] = elStiff_i[jcol++];
84 copy[j][i] = copy[i][j];
85 }
86 }
87 break;
88
89 case FEI_LOWER_SYMM_ROW:
90 for (i = 0; i < numRows; i++) {
91 elStiff_i = elemStiff[i];
92 for (j = 0; j <=i; j++) {
93 copy[i][j] = elStiff_i[j];
94 copy[j][i] = copy[i][j];
95 }
96 }
97 break;
98
99 case FEI_DENSE_COL:
100 for (i = 0; i < numRows; i++) {
101 elStiff_i = elemStiff[i];
102 for (j = 0; j < numRows; j++) {
103 copy[j][i] = elStiff_i[j];
104 }
105 }
106 break;
107
108 case FEI_UPPER_SYMM_COL:
109 for (i = 0; i < numRows; i++) {
110 elStiff_i = elemStiff[i];
111 for (j = 0; j <= i; j++) {
112 copy[i][j] = elStiff_i[j];
113 copy[j][i] = copy[i][j];
114 }
115 }
116 break;
117
118 case FEI_LOWER_SYMM_COL:
119 for (i = 0; i < numRows; i++) {
120 elStiff_i = elemStiff[i];
121 int jcol=0;
122 for (j = i; j < numRows; j++) {
123 copy[i][j] = elStiff_i[jcol++];
124 copy[j][i] = copy[i][j];
125 }
126 }
127 break;
128
129 default:
130 throw std::runtime_error("copyStiffness ERROR, unrecognized elem-format");
131 }
132}
133
134//------------------------------------------------------------------------------
135const NodeDescriptor* Filter::findNode(GlobalID nodeID) const {
136//
137//This function returns a NodeDescriptor ptr, may return NULL.
138//
139 const NodeDescriptor* node = NULL;
140 problemStructure_->getNodeDatabase().getNodeWithID(nodeID, node);
141 return node;
142}
143
144//------------------------------------------------------------------------------
145const NodeDescriptor& Filter::findNodeDescriptor(GlobalID nodeID) const {
146//
147//This function returns a NodeDescriptor reference if nodeID is an active node.
148//
149 const NodeDescriptor* node = NULL;
150 int err = problemStructure_->getNodeDatabase().getNodeWithID(nodeID, node);
151
152 if (err != 0) {
153 fei::console_out() << "ERROR, Filter::findNodeDescriptor unable to find node "
154 << static_cast<int>(nodeID) << FEI_ENDL;
155 std::abort();
156 }
157
158 return( *node );
159}
160//------------------------------------------------------------------------------
161int Filter::calculateResidualNorms(int whichNorm, int numFields,
162 int* fieldIDs, double* norms,
163 std::vector<double>& residValues)
164{
165 std::vector<double> normsArray(numFields, 0.0);
166
167 std::fill(fieldIDs, fieldIDs+numFields, -999);
168
169 std::vector<double> tmpNorms(numFields);
170 double* tmpNormsPtr = &tmpNorms[0];
171
172 double* residPtr = &(residValues[0]);
173
174 const std::vector<int>& pfieldIDs = problemStructure_->getFieldIDs();
175 int numDBFields = pfieldIDs.size();
176 std::vector<int>::const_iterator
177 f_iter = pfieldIDs.begin(),
178 f_end = pfieldIDs.end();
179
180 int DBFieldSize = 0;
181
182 int offset = 0;
183 for(; f_iter != f_end; ++f_iter) {
184
185 if (offset == 0) DBFieldSize = problemStructure_->getFieldSize(*f_iter);
186
187 if (*f_iter > -1) {
188 if (offset < numFields) {
189 fieldIDs[offset] = *f_iter;
190 tmpNormsPtr[offset++] = 0.0;
191 }
192 }
193 }
194
195 int reducedStartRow = problemStructure_->getFirstReducedEqn();
196 int reducedEndRow = problemStructure_->getLastReducedEqn();
197
198 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
199 int numNodes = nodeDB.getNumNodeDescriptors();
200
201 bool haveSlaves = problemStructure_->numSlaveEquations() > 0;
202
203 for(int i=0; i<numNodes; i++) {
204 const NodeDescriptor* node = NULL;
205 nodeDB.getNodeAtIndex(i, node);
206
207 if (node==NULL || node->getOwnerProc() != localRank_) continue;
208
209 const int* fieldIDList = node->getFieldIDList();
210 const int* fieldEqnNums = node->getFieldEqnNumbers();
211 int numNodeFields = node->getNumFields();
212
213 for(int j=0; j<numNodeFields; j++) {
214 int fIndex = 0;
215 int fSize = DBFieldSize;
216
217 if (numDBFields > 1) {
218 fIndex = fei::binarySearch(fieldIDList[j], fieldIDs, numFields);
219 if (fIndex < 0) return(-1);
220 fSize = problemStructure_->getFieldSize(fieldIDList[j]);
221 if (fSize < 0) return(-1);
222 }
223
224 for(int k=0; k<fSize; k++) {
225 int eqn = fieldEqnNums[j]+k;
226
227 if (haveSlaves) {
228 if (problemStructure_->isSlaveEqn(eqn)) continue;
229 int reducedEqn;
230 problemStructure_->translateToReducedEqn(eqn, reducedEqn);
231
232 if (reducedEqn < reducedStartRow || reducedEqn > reducedEndRow) {
233 continue;
234 }
235 eqn = reducedEqn;
236 }
237
238 double rval = residPtr[eqn - reducedStartRow];
239
240 switch(whichNorm) {
241 case 0:
242 if (tmpNormsPtr[fIndex] < std::abs(rval)) tmpNormsPtr[fIndex] = rval;
243 break;
244 case 1:
245 tmpNormsPtr[fIndex] += std::abs(rval);
246 break;
247 case 2:
248 tmpNormsPtr[fIndex] += rval*rval;
249 break;
250 default:
251 FEI_COUT << "Filter::residualNorm: ERROR, whichNorm="<<whichNorm
252 << " not recognized." << FEI_ENDL;
253 return(-1);
254 }
255 }
256 }
257 }
258
259 //so at this point we have the local norms. We now need to perform the
260 //global max or sum, depending on which norm it is.
261
262 MPI_Comm comm = problemStructure_->getCommunicator();
263
264 if (whichNorm != 0) {
265 CHK_ERR( fei::GlobalSum(comm, tmpNorms, normsArray) );
266 }
267 else {
268 CHK_ERR( fei::GlobalMax(comm, tmpNorms, normsArray) );
269 }
270
271 for(int i=0; i<numFields; ++i) {
272 norms[i] = normsArray[i];
273 }
274
275 if (whichNorm == 2) {
276 for(int i=0; i<numFields; ++i) norms[i] = std::sqrt(norms[i]);
277 }
278
279 return(0);
280}
281
282//------------------------------------------------------------------------------
283int Filter::parameters(int numParams, const char *const* paramStrings)
284{
285 if (numParams == 0 || paramStrings == NULL) return(0);
286
287 const char* param = snl_fei::getParam("outputLevel",numParams,paramStrings);
288
289 if ( param != NULL){
290 std::string str(&(param[11]));
291 FEI_ISTRINGSTREAM isstr(str);
292 isstr >> outputLevel_;
293 }
294
295 return(0);
296}
297
Filter(SNL_FEI_Structure *probStruct)
virtual ~Filter()
int getNodeWithID(GlobalID nodeID, const NodeDescriptor *&node) const
void getNodeAtIndex(int i, const NodeDescriptor *&node) const
int getNumNodeDescriptors() const
int getFieldSize(int fieldID)
bool translateToReducedEqn(int eqn, int &reducedEqn)
int binarySearch(const T &item, const T *list, int len)
std::ostream & console_out()
int GlobalMax(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
int GlobalSum(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)