FEI Version of the Day
Loading...
Searching...
No Matches
fei_LinSysCoreFilter.hpp
1#ifndef _LinSysCoreFilter_hpp_
2#define _LinSysCoreFilter_hpp_
3
4/*--------------------------------------------------------------------*/
5/* Copyright 2005 Sandia Corporation. */
6/* Under the terms of Contract DE-AC04-94AL85000, there is a */
7/* non-exclusive license for use of this work by or on behalf */
8/* of the U.S. Government. Export of this program may require */
9/* a license from the United States Government. */
10/*--------------------------------------------------------------------*/
11
12#include "fei_macros.hpp"
13#include "fei_defs.h"
14#include "fei_fwd.hpp"
15#include "fei_Filter.hpp"
16#include <fei_CSRMat.hpp>
17#include <fei_CSVec.hpp>
18
19namespace fei {
20 class DirichletBCManager;
21}
22
33class LinSysCoreFilter : public Filter {
34
35 public:
36 // Constructor.
37 LinSysCoreFilter(FEI_Implementation* owner, MPI_Comm comm,
38 SNL_FEI_Structure* probStruct,
40 int masterRank=0);
41
42 //Destructor
43 virtual ~LinSysCoreFilter();
44
45
46 // set a value (usually zeros) throughout the linear system
47 virtual int resetSystem(double s);
48 virtual int resetMatrix(double s);
49 virtual int resetRHSVector(double s);
50 virtual int resetInitialGuess(double s);
51
52 virtual int deleteMultCRs();
53
54 virtual int loadNodeBCs(int numNodes,
55 const GlobalID *nodeIDs,
56 int fieldID,
57 const int* offsetsIntoField,
58 const double* prescribedValues);
59
60 virtual int loadElemBCs(int numElems,
61 const GlobalID *elemIDs,
62 int fieldID,
63 const double *const *alpha,
64 const double *const *beta,
65 const double *const *gamma);
66
67 virtual int sumInElem(GlobalID elemBlockID,
68 GlobalID elemID,
69 const GlobalID* elemConn,
70 const double* const* elemStiffness,
71 const double* elemLoad,
72 int elemFormat);
73
74 virtual int sumInElemMatrix(GlobalID elemBlockID,
75 GlobalID elemID,
76 const GlobalID* elemConn,
77 const double* const* elemStiffness,
78 int elemFormat);
79
80 virtual int sumInElemRHS(GlobalID elemBlockID,
81 GlobalID elemID,
82 const GlobalID* elemConn,
83 const double* elemLoad);
84
85 virtual int loadCRMult(int CRMultID,
86 int numCRNodes,
87 const GlobalID* CRNodes,
88 const int* CRFields,
89 const double* CRWeights,
90 double CRValue);
91
92 virtual int loadCRPen(int CRPenID,
93 int numCRNodes,
94 const GlobalID* CRNodes,
95 const int *CRFields,
96 const double* CRWeights,
97 double CRValue,
98 double penValue);
99
100 virtual int putIntoRHS(int IDType,
101 int fieldID,
102 int numIDs,
103 const GlobalID* IDs,
104 const double* rhsEntries);
105
106 virtual int sumIntoRHS(int IDType,
107 int fieldID,
108 int numIDs,
109 const GlobalID* IDs,
110 const double* rhsEntries);
111
112 virtual int loadComplete();
113
114 // set parameters associated with solver choice, etc.
115 virtual int parameters(int numParams, const char *const* paramStrings);
116
117 //get residual norms
118 virtual int residualNorm(int whichNorm, int numFields,
119 int* fieldIDs, double* norms, double& residTime);
120
121 // start iterative solution
122 virtual int solve(int& status, double& sTime);
123
124 // query function iterations performed.
125 virtual int iterations() const {return(iterations_);};
126
127// Solution return services.......................................
128
129 // return all nodal solution params on a block-by-block basis
130 virtual int getBlockNodeSolution(GlobalID elemBlockID,
131 int numNodes,
132 const GlobalID *nodeIDs,
133 int *offsets,
134 double *results);
135
136 virtual int getNodalSolution(int numNodes,
137 const GlobalID *nodeIDs,
138 int *offsets,
139 double *results);
140
141 // return nodal solution for one field on a block-by-block basis
142 virtual int getBlockFieldNodeSolution(GlobalID elemBlockID,
143 int fieldID,
144 int numNodes,
145 const GlobalID *nodeIDs,
146 double *results);
147
148 // return element solution params on a block-by-block basis
149 virtual int getBlockElemSolution(GlobalID elemBlockID,
150 int numElems,
151 const GlobalID *elemIDs,
152 int& numElemDOFPerElement,
153 double *results);
154
155 virtual int getCRMultipliers(int numCRs, const int* CRIDs, double* multipliers);
156
157// associated "puts" paralleling the solution return services.
158//
159// the int sizing parameters are passed for error-checking purposes, so
160// that the interface implementation can tell if the passed estimate
161// vectors make sense -before- an attempt is made to utilize them as
162// initial guesses by unpacking them into the solver's native solution
163// vector format (these parameters include lenNodeIDList, lenElemIDList,
164// numElemDOF, and numMultCRs -- all other passed params are either
165// vectors or block/constraint-set IDs)
166
167 // put nodal-based solution guess on a block-by-block basis
168 virtual int putBlockNodeSolution(GlobalID elemBlockID,
169 int numNodes,
170 const GlobalID *nodeIDs,
171 const int *offsets,
172 const double *estimates);
173
174 // put nodal-based guess for one field on a block-by-block basis
175 virtual int putBlockFieldNodeSolution(GlobalID elemBlockID,
176 int fieldID,
177 int numNodes,
178 const GlobalID *nodeIDs,
179 const double *estimates);
180
181 // put element-based solution guess on a block-by-block basis
182 virtual int putBlockElemSolution(GlobalID elemBlockID,
183 int numElems,
184 const GlobalID *elemIDs,
185 int dofPerElem,
186 const double *estimates);
187
188 virtual int putCRMultipliers(int numMultCRs,
189 const int* CRIDs,
190 const double *multEstimates);
191
192//===== a couple of public non-FEI functions... ================================
193//These are intended to be used by an 'outer-layer' class like
194//FEI_Implementation.
195//
196 public:
197 virtual int getNodalFieldSolution(int fieldID,
198 int numNodes,
199 const GlobalID* nodeIDs,
200 double* results);
201
202 virtual int putNodalFieldData(int fieldID,
203 int numNodes,
204 const GlobalID* nodeIDs,
205 const double* nodeData);
206
207 virtual int putNodalFieldSolution(int fieldID,
208 int numNodes,
209 const GlobalID* nodeIDs,
210 const double* nodeData);
211
212 virtual int unpackSolution();
213
214 void setEqnCommMgr(EqnCommMgr* eqnCommMgr);
215
216 EqnCommMgr* getEqnCommMgr() {return(eqnCommMgr_);};
217
218 virtual int setNumRHSVectors(int numRHSs, int* rhsIDs);
219 virtual int setCurrentRHS(int rhsID);
220
221 virtual int exchangeRemoteEquations();
222 virtual int exchangeRemoteBCs(std::vector<int>& essEqns,
223 std::vector<double>& essAlpha,
224 std::vector<double>& essGamma);
225
226 virtual int implementAllBCs();
227
228 virtual int enforceEssentialBCs(const int* eqns, const double* alpha,
229 const double* gamma, int numEqns);
230
231 virtual int enforceRemoteEssBCs(int numEqns, const int* eqns,
232 const int* const* colIndices, const int* colIndLens,
233 const double* const* coefs);
234
235 virtual int initialize();
236
237//==============================================================================
238//private functions for internal implementation of LinSysCoreFilter.
239//==============================================================================
240 private:
241 int initLinSysCore();
242 void setLinSysCoreCREqns();
243
244 int unpackRemoteContributions(EqnCommMgr& eqnCommMgr,
245 int assemblyMode);
246
247 int loadFEDataMultCR(int CRID,
248 int numCRNodes,
249 const GlobalID* CRNodes,
250 const int* CRFields,
251 const double* CRWeights,
252 double CRValue);
253
254 int loadFEDataPenCR(int CRID,
255 int numCRNodes,
256 const GlobalID* CRNodes,
257 const int* CRFields,
258 const double* CRWeights,
259 double CRValue,
260 double penValue);
261
262 int storeNodalColumnCoefs(int eqn, const NodeDescriptor& node,
263 int fieldID, int fieldSize,
264 double* coefs);
265
266 int storeNodalRowCoefs(const NodeDescriptor& node,
267 int fieldID, int fieldSize,
268 double* coefs, int eqn);
269
270 int generalElemInput(GlobalID elemBlockID,
271 GlobalID elemID,
272 const double* const* elemStiffness,
273 const double* elemLoad,
274 int elemFormat);
275
276 int generalElemInput(GlobalID elemBlockID,
277 GlobalID elemID,
278 const GlobalID* elemConn,
279 const double* const* elemStiffness,
280 const double* elemLoad,
281 int elemFormat);
282
283 void storeNodalSendIndex(const NodeDescriptor& node, int fieldID, int col);
284 void storeNodalSendEqn(const NodeDescriptor& node, int fieldID, int col,
285 double* coefs);
286 void storeNodalSendIndices(const NodeDescriptor& iNode, int iField,
287 const NodeDescriptor& jNode, int jField);
288
289 void storePenNodeSendData(const NodeDescriptor& iNode,
290 int iField, int iFieldSize,
291 double* iCoefs,
292 const NodeDescriptor& jNode,
293 int jField, int jFieldSize,
294 double* jCoefs,
295 double penValue, double CRValue);
296
297 int storePenNodeData(const NodeDescriptor& iNode,
298 int iField, int iFieldSize,
299 double* iCoefs,
300 const NodeDescriptor& jNode,
301 int jField, int jFieldSize,
302 double* jCoefs,
303 double penValue, double CRValue);
304
305 void allocElemStuff();
306
307 int resolveConflictingCRs(EqnBuffer& bcEqns);
308
309 int giveToMatrix_symm_noSlaves(int numPtRows,
310 const int* ptRowNumbers,
311 const double* const* coefs,
312 int mode);
313
314 int giveToBlkMatrix_symm_noSlaves(int numPtRows, const int* ptRows,
315 int numBlkRows, const int* blkRowNumbers,
316 const int* blkRowSizes,
317 const double* const* coefs,
318 int mode);
319
320 int giveToMatrix(int numPtRows, const int* ptRows,
321 int numPtCols, const int* ptCols,
322 const double* const* values,
323 int mode);
324
325 int giveToLocalReducedMatrix(int numPtRows, const int* ptRows,
326 int numPtCols, const int* ptCols,
327 const double* const* values,
328 int mode);
329
330 int getFromMatrix(int numPtRows, const int* ptRows,
331 const int* rowColOffsets, const int* ptCols,
332 int numColsPerRow, double** values);
333
334 int sumIntoMatrix(fei::CSRMat& mat);
335
336 int getEqnsFromMatrix(ProcEqns& procEqns, EqnBuffer& eqnData);
337
338 int getEqnsFromRHS(ProcEqns& procEqns, EqnBuffer& eqnData);
339
340 int giveToRHS(int num, const double* values,
341 const int* indices, int mode);
342
343 int giveToLocalReducedRHS(int num, const double* values,
344 const int* indices, int mode);
345
346 int getFromRHS(int num, double* values, const int* indices);
347
348 int sumIntoRHS(fei::CSVec& vec);
349
350 int getEqnSolnEntry(int eqnNumber, double& solnValue);
351
352 int getSharedRemoteSolnEntry(int eqnNumber, double& solnValue);
353
354 int getReducedSolnEntry(int eqnNumber, double& solnValue);
355
356 int formResidual(double* residValues, int numLocalEqns);
357
358 int getRemoteSharedEqns(int numPtRows, const int* ptRows,
359 ProcEqns& remoteProcEqns);
360
361 int resetTheMatrix(double s);
362 int resetTheRHSVector(double s);
363
364 int assembleEqns(int numPtRows,
365 int numPtCols,
366 const int* rowNumbers,
367 const int* colIndices,
368 const double* const* coefs,
369 bool structurallySymmetric,
370 int numBlkEqns, int* blkEqns, int* blkSizes,
371 bool useBlkEqns, int mode);
372
373 int assembleReducedEqns();
374
375 int assembleRHS(int numValues, const int* indices, const double* coefs, int mode);
376
377 int assembleReducedRHS();
378
379 void debugOutput(const char* mesg);
380
381 int createEqnCommMgr_put();
382
383//==============================================================================
384//private LinSysCoreFilter variables
385//==============================================================================
386 private:
387
388 int timesInitializeCalled_;
389
390 LinearSystemCore* lsc_;
391 bool useLookup_;
392
393 int internalFei_;
394
395 bool newMatrixData_;
396 bool newVectorData_;
397 bool newConstraintData_;
398 bool newBCData_;
399 bool connectivitiesInitialized_;
400 bool firstRemEqnExchange_;
401 bool needToCallMatrixLoadComplete_;
402
403 bool resolveConflictRequested_;
404
405 int localStartRow_, localEndRow_, numLocalEqns_, numGlobalEqns_;
406 int reducedStartRow_, reducedEndRow_, numReducedRows_;
407 int numLocallyOwnedNodes_, numGlobalNodes_, firstLocalNodeNumber_;
408
409 bool blockMatrix_;
410 bool tooLateToChooseBlock_;
411
412 int numLocalEqnBlks_;
413 int localReducedBlkOffset_, numLocalReducedEqnBlks_;
414
415 int iterations_;
416 int numRHSs_;
417 int currentRHS_;
418 std::vector<int> rhsIDs_;
419
420 int outputLevel_;
421
422 MPI_Comm comm_;
423 int masterRank_;
424
425 SNL_FEI_Structure* problemStructure_;
426 bool matrixAllocated_;
427
428 std::vector<int> rowIndices_;
429 std::vector<int> rowColOffsets_, colIndices_;
430
431 fei::FillableMat *Kid_, *Kdi_, *Kdd_;
432 fei::CSRMat csrD, csrKid, csrKdi, csrKdd, tmpMat1_, tmpMat2_;
433 fei::CSVec fd_, tmpVec1_;
434 int reducedEqnCounter_, reducedRHSCounter_;
435 std::vector<int> rSlave_, cSlave_;
436
437 int nodeIDType_;
438 fei::DirichletBCManager* bcManager_; //Boundary condition manager
439
440 EqnCommMgr* eqnCommMgr_; //equation communication manager
441 EqnCommMgr* eqnCommMgr_put_; //only created if users call the 'put'
442 // functions
443
444 int maxElemRows_;
445 std::vector<int> scatterIndices_;
446 std::vector<int> blkScatterIndices_;
447 std::vector<int> iworkSpace_, iworkSpace2_;
448 std::vector<double> dworkSpace_;
449 std::vector<const double*> dworkSpace2_;
450
451 double** eStiff_;
452 double* eStiff1D_;
453 double* eLoad_;
454};
455
456#endif
457