FEI Version of the Day
Loading...
Searching...
No Matches
snl_fei_LinearSystem_General.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#include <fei_utils.hpp>
11
12#include <cmath>
13
14#include <snl_fei_LinearSystem_General.hpp>
15#include <fei_MatrixReducer.hpp>
16#include <fei_Matrix_Impl.hpp>
17#include <fei_VectorSpace.hpp>
18#include <fei_MatrixGraph.hpp>
19#include <fei_SparseRowGraph.hpp>
20#include <snl_fei_Constraint.hpp>
21#include <fei_Record.hpp>
22#include <fei_utils.hpp>
23#include <fei_impl_utils.hpp>
24#include <fei_LogManager.hpp>
25
26#include <fei_DirichletBCRecord.hpp>
27#include <fei_DirichletBCManager.hpp>
28#include <fei_EqnBuffer.hpp>
29#include <fei_LinSysCoreFilter.hpp>
30
31#undef fei_file
32#define fei_file "snl_fei_LinearSystem_General.cpp"
33#include <fei_ErrMacros.hpp>
34
35//----------------------------------------------------------------------------
37 : fei::LinearSystem(matrixGraph),
38 comm_(matrixGraph->getRowSpace()->getCommunicator()),
39 essBCvalues_(NULL),
40 resolveConflictRequested_(false),
41 bcs_trump_slaves_(false),
42 explicitBCenforcement_(false),
43 BCenforcement_no_column_mod_(false),
44 localProc_(0),
45 numProcs_(1),
46 name_(),
47 named_loadcomplete_counter_(),
48 iwork_(),
49 dwork_(),
50 dbgprefix_("LinSysG: ")
51{
52 localProc_ = fei::localProc(comm_);
53 numProcs_ = fei::numProcs(comm_);
54
55 fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph->getRowSpace();
56
57 std::vector<int> offsets;
58 vecSpace->getGlobalIndexOffsets(offsets);
59
60 firstLocalOffset_ = offsets[localProc_];
61 lastLocalOffset_ = offsets[localProc_+1]-1;
62
63 setName("dbg");
64}
65
66//----------------------------------------------------------------------------
68{
69 delete essBCvalues_;
70}
71
72//----------------------------------------------------------------------------
74 const char* const* paramStrings)
75{
76 if (numParams == 0 || paramStrings == NULL) return(0);
77
78 const char* param = snl_fei::getParam("name", numParams, paramStrings);
79 if (param != NULL) {
80 if (strlen(param) < 6) ERReturn(-1);
81
82 setName(&(param[5]));
83 }
84
85 param = snl_fei::getParam("resolveConflict",numParams,paramStrings);
86 if (param != NULL){
87 resolveConflictRequested_ = true;
88 }
89
90 param = snl_fei::getParam("BCS_TRUMP_SLAVE_CONSTRAINTS",
91 numParams,paramStrings);
92 if (param != NULL) {
93 bcs_trump_slaves_ = true;
94 }
95
96 param = snl_fei::getParam("EXPLICIT_BC_ENFORCEMENT",numParams,paramStrings);
97 if (param != NULL){
98 explicitBCenforcement_ = true;
99 }
100
101 param = snl_fei::getParam("BC_ENFORCEMENT_NO_COLUMN_MOD",numParams,paramStrings);
102 if (param != NULL){
103 BCenforcement_no_column_mod_ = true;
104 }
105
106 param = snl_fei::getParamValue("FEI_OUTPUT_LEVEL",numParams,paramStrings);
107 if (param != NULL) {
108 setOutputLevel(fei::utils::string_to_output_level(param));
109 }
110
111 if (matrix_.get() != NULL) {
112 fei::Matrix* matptr = matrix_.get();
113 fei::MatrixReducer* matred = dynamic_cast<fei::MatrixReducer*>(matptr);
114 if (matred != NULL) {
115 matptr = matred->getTargetMatrix().get();
116 }
118 dynamic_cast<fei::Matrix_Impl<LinearSystemCore>*>(matptr);
119 if (lscmatrix != NULL) {
120 lscmatrix->getMatrix()->parameters(numParams, (char**)paramStrings);
121 }
122 }
123
124 return(0);
125}
126
127//----------------------------------------------------------------------------
129{
130 int numParams = 0;
131 const char** paramStrings = NULL;
132 std::vector<std::string> stdstrings;
134 fei::utils::strings_to_char_ptrs(stdstrings, numParams, paramStrings);
135
136 int err = parameters(numParams, paramStrings);
137
138 delete [] paramStrings;
139
140 return(err);
141}
142
143//----------------------------------------------------------------------------
144void snl_fei::LinearSystem_General::setName(const char* name)
145{
146 if (name == NULL) return;
147
148 if (name_ == name) return;
149
150 name_ = name;
151
152 std::map<std::string,unsigned>::iterator
153 iter = named_loadcomplete_counter_.find(name_);
154 if (iter == named_loadcomplete_counter_.end()) {
155 static int counter = 0;
156 named_loadcomplete_counter_.insert(std::make_pair(name_, counter));
157 ++counter;
158 }
159}
160
161//----------------------------------------------------------------------------
163 const int* IDs,
164 int idType,
165 int fieldID,
166 int offsetIntoField,
167 const double* prescribedValues)
168{
169 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != 0) {
170 FEI_OSTREAM& os = *output_stream_;
171 os << "loadEssentialBCs, numIDs: "<<numIDs<<", idType: " <<idType
172 <<", fieldID: "<<fieldID<<", offsetIntoField: "<<offsetIntoField<<FEI_ENDL;
173 }
174
175 return fei::LinearSystem::loadEssentialBCs(numIDs, IDs, idType, fieldID,
176 offsetIntoField, prescribedValues);
177}
178
179//----------------------------------------------------------------------------
181 const int* IDs,
182 int idType,
183 int fieldID,
184 const int* offsetsIntoField,
185 const double* prescribedValues)
186{
187 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != 0) {
188 FEI_OSTREAM& os = *output_stream_;
189 for(int i=0; i<numIDs; ++i) {
190 os << "loadEssentialBCs idType: " <<idType
191 <<", fieldID: "<<fieldID<<", ID: " << IDs[i]<<", offsetIntoField: "<<offsetsIntoField[i]<<", val: " << prescribedValues[i] << FEI_ENDL;
192 }
193 }
194
195 return fei::LinearSystem::loadEssentialBCs(numIDs, IDs, idType, fieldID,
196 offsetsIntoField, prescribedValues);
197}
198
199//----------------------------------------------------------------------------
201 bool globalAssemble)
202{
203 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != 0) {
204 FEI_OSTREAM& os = *output_stream_;
205 os << dbgprefix_<<"loadComplete"<<FEI_ENDL;
206 }
207
208 if (dbcManager_ == NULL) {
209 dbcManager_ = new fei::DirichletBCManager(matrixGraph_->getRowSpace());
210 }
211
212 if (globalAssemble) {
213
214 if (matrix_.get() != NULL) {
215 CHK_ERR( matrix_->gatherFromOverlap() );
216 }
217
218 if (rhs_.get() != NULL) {
219 CHK_ERR( rhs_->gatherFromOverlap() );
220 }
221
222 }
223
224 unsigned counter = 0;
225
226 std::map<std::string,unsigned>::iterator
227 iter = named_loadcomplete_counter_.find(name_);
228 if (iter == named_loadcomplete_counter_.end()) {
229 FEI_COUT << "fei::LinearSystem::loadComplete internal error, name "
230 << name_ << " not found." << FEI_ENDL;
231 }
232 else {
233 counter = iter->second++;
234 }
235
236 if (output_level_ >= fei::FULL_LOGS) {
237 std::string opath = fei::LogManager::getLogManager().getOutputPath();
238 if (opath == "") opath = ".";
239
240 FEI_OSTRINGSTREAM Aname;
241 FEI_OSTRINGSTREAM bname;
242 FEI_OSTRINGSTREAM xname;
243 Aname << opath << "/";
244 bname << opath << "/";
245 xname << opath << "/";
246
247 Aname << "A_"<<name_<<".preBC.np"<<numProcs_<<".slv"<<counter<< ".mtx";
248
249 bname << "b_"<<name_<<".preBC.np"<<numProcs_<<".slv"<<counter<< ".vec";
250
251 std::string Aname_str = Aname.str();
252 const char* Aname_c_str = Aname_str.c_str();
253 CHK_ERR( matrix_->writeToFile(Aname_c_str) );
254
255 std::string bname_str = bname.str();
256 const char* bname_c_str = bname_str.c_str();
257 CHK_ERR( rhs_->writeToFile(bname_c_str) );
258 }
259
260 CHK_ERR( implementBCs(applyBCs) );
261
262 if (globalAssemble) {
263 CHK_ERR( matrix_->globalAssemble() );
264 }
265
266 if (output_level_ == fei::STATS || output_level_ == fei::ALL) {
267 int globalNumSlaveCRs = matrixGraph_->getGlobalNumSlaveConstraints();
268 if (localProc_ == 0) {
269 FEI_COUT << "Global Neqns: " << matrix_->getGlobalNumRows();
270 if (globalNumSlaveCRs > 0) {
271 FEI_COUT << ", Global NslaveCRs: " << globalNumSlaveCRs;
272 }
273 FEI_COUT << FEI_ENDL;
274 }
275 }
276
277 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
278 FEI_OSTREAM& os = *output_stream_;
279 os << dbgprefix_<<"Neqns=" << matrix_->getGlobalNumRows();
280 int globalNumSlaveCRs = matrixGraph_->getGlobalNumSlaveConstraints();
281 if (globalNumSlaveCRs > 0) {
282 os << ", Global NslaveCRs=" << globalNumSlaveCRs;
283 }
284 os << FEI_ENDL;
285 }
286
287 if (output_level_ >= fei::MATRIX_FILES) {
288 std::string opath = fei::LogManager::getLogManager().getOutputPath();
289 if (opath == "") opath = ".";
290
291 FEI_OSTRINGSTREAM Aname;
292 FEI_OSTRINGSTREAM bname;
293 FEI_OSTRINGSTREAM xname;
294 Aname << opath << "/";
295 bname << opath << "/";
296 xname << opath << "/";
297
298 Aname << "A_" <<name_<<".np"<<numProcs_<< ".slv" << counter << ".mtx";
299
300 bname << "b_" <<name_<<".np"<<numProcs_<< ".slv" << counter << ".vec";
301
302 xname << "x0_" <<name_<<".np"<<numProcs_<< ".slv" << counter << ".vec";
303
304 std::string Aname_str = Aname.str();
305 const char* Aname_c_str = Aname_str.c_str();
306 CHK_ERR( matrix_->writeToFile(Aname_c_str) );
307
308 std::string bname_str = bname.str();
309 const char* bname_c_str = bname_str.c_str();
310 CHK_ERR( rhs_->writeToFile(bname_c_str) );
311
312 std::string xname_str = xname.str();
313 const char* xname_c_str = xname_str.c_str();
314 CHK_ERR( soln_->writeToFile(xname_c_str) );
315 }
316
317 return(0);
318}
319
320//----------------------------------------------------------------------------
322{
323 if (essBCvalues_ == NULL) {
324 return(0);
325 }
326
327 if (essBCvalues_->size() == 0) {
328 return(0);
329 }
330
331 CHK_ERR( vector->copyIn(essBCvalues_->size(),
332 &(essBCvalues_->indices())[0],
333 &(essBCvalues_->coefs())[0]) );
334
335 return(0);
336}
337
338//----------------------------------------------------------------------------
340{
341 if (essBCvalues_ == NULL) return(false);
342
343 std::vector<int>& indices = essBCvalues_->indices();
344 int offset = fei::binarySearch(globalEqnIndex, indices);
345 return( offset < 0 ? false : true);
346}
347
348//----------------------------------------------------------------------------
350 std::vector<double>& bcVals) const
351{
352 bcEqns.clear();
353 bcVals.clear();
354 if (essBCvalues_ == NULL) return;
355
356 int num = essBCvalues_->size();
357 bcEqns.resize(num);
358 bcVals.resize(num);
359 int* essbcs = &(essBCvalues_->indices())[0];
360 double* vals = &(essBCvalues_->coefs())[0];
361 for(int i=0; i<num; ++i) {
362 bcEqns[i] = essbcs[i];
363 bcVals[i] = vals[i];
364 }
365}
366
367//----------------------------------------------------------------------------
368void snl_fei::LinearSystem_General::getConstrainedEqns(std::vector<int>& crEqns) const
369{
370 matrixGraph_->getConstrainedIndices(crEqns);
371}
372
373//----------------------------------------------------------------------------
374int extractDirichletBCs(fei::DirichletBCManager* bcManager,
376 fei::CSVec* essBCvalues,
377 bool resolveConflictRequested,
378 bool bcs_trump_slaves)
379{
380// int numLocalBCs = bcManager->getNumBCRecords();
381// int globalNumBCs = 0;
382// MPI_Comm comm = matrixGraph->getRowSpace()->getCommunicator();
383// fei::GlobalSum(comm, numLocalBCs, globalNumBCs);
384// if (globalNumBCs == 0) {
385// return(0);
386// }
387
388 fei::SharedPtr<fei::FillableMat> localBCeqns(new fei::FillableMat);
390// matrixGraph->getRowSpace()->initComplete();
391 int numSlaves = matrixGraph->getGlobalNumSlaveConstraints();
392 fei::SharedPtr<fei::Reducer> reducer = matrixGraph->getReducer();
393
394 int numIndices = numSlaves>0 ?
395 reducer->getLocalReducedEqns().size() :
396 matrixGraph->getRowSpace()->getNumIndices_Owned();
397
398 bool zeroSharedRows = false;
399 bcEqns.reset(new fei::Matrix_Impl<fei::FillableMat>(localBCeqns, matrixGraph, numIndices, zeroSharedRows));
400 fei::SharedPtr<fei::Matrix> bcEqns_reducer;
401 if (numSlaves > 0) {
402 bcEqns_reducer.reset(new fei::MatrixReducer(reducer, bcEqns));
403 }
404
405 fei::Matrix& bcEqns_mat = bcEqns_reducer.get()==NULL ?
406 *bcEqns : *bcEqns_reducer;
407
408 CHK_ERR( bcManager->finalizeBCEqns(bcEqns_mat, bcs_trump_slaves) );
409
410 if (resolveConflictRequested) {
411 fei::SharedPtr<fei::FillableMat> mat = bcEqns->getMatrix();
412 std::vector<int> bcEqnNumbers;
413 fei::get_row_numbers(*mat, bcEqnNumbers);
414 CHK_ERR( snl_fei::resolveConflictingCRs(*matrixGraph, bcEqns_mat,
415 bcEqnNumbers) );
416 }
417
418 std::vector<int> essEqns;
419 std::vector<double> values;
420
421 std::map<int,fei::FillableMat*>& remotes = bcEqns->getRemotelyOwnedMatrices();
422 std::map<int,fei::FillableMat*>::iterator
423 it = remotes.begin(),
424 it_end = remotes.end();
425 for(; it!=it_end; ++it) {
426 fei::impl_utils::separate_BC_eqns( *(it->second), essEqns, values);
427 }
428
429// CHK_ERR( bcEqns->gatherFromOverlap(false) );
430
431 fei::impl_utils::separate_BC_eqns( *(bcEqns->getMatrix()), essEqns, values);
432
433 if (essEqns.size() > 0) {
434 int* essEqnsPtr = &essEqns[0];
435 double* valuesPtr = &values[0];
436
437 for(unsigned i=0; i<essEqns.size(); ++i) {
438 int eqn = essEqnsPtr[i];
439 double value = valuesPtr[i];
440 fei::put_entry(*essBCvalues, eqn, value);
441 }
442 }
443
444 return(0);
445}
446
447//----------------------------------------------------------------------------
448int snl_fei::LinearSystem_General::implementBCs(bool applyBCs)
449{
450 if (essBCvalues_ != NULL) {
451 delete essBCvalues_;
452 }
453
454 essBCvalues_ = new fei::CSVec;
455
456 CHK_ERR( extractDirichletBCs(dbcManager_, matrixGraph_,
457 essBCvalues_, resolveConflictRequested_,
458 bcs_trump_slaves_) );
459
460 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
461 FEI_OSTREAM& os = *output_stream_;
462 std::vector<int>& indices = essBCvalues_->indices();
463 std::vector<double>& coefs= essBCvalues_->coefs();
464 for(size_t i=0; i<essBCvalues_->size(); ++i) {
465 os << "essBCeqns["<<i<<"]: "<<indices[i]<<", "<<coefs[i]<<FEI_ENDL;
466 }
467 }
468
469 //If the underlying matrix is a LinearSystemCore instance, then this
470 //function will return 0, and we're done. A non-zero return-code means
471 //we should continue and enforce the BCs assuming a general matrix.
472
473 int returncode = enforceEssentialBC_LinSysCore();
474 if (returncode == 0) {
475 return(0);
476 }
477
478 fei::CSVec allEssBCs;
479 if (!BCenforcement_no_column_mod_) {
480 fei::impl_utils::global_union(comm_, *essBCvalues_, allEssBCs);
481
482 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
483 FEI_OSTREAM& os = *output_stream_;
484 os << " implementBCs, essBCvalues_.size(): "<<essBCvalues_->size()
485 << ", allEssBCs.size(): " << allEssBCs.size()<<FEI_ENDL;
486 }
487 }
488
489 if (essBCvalues_->size() > 0) {
490 enforceEssentialBC_step_1(*essBCvalues_);
491 }
492
493 if (!BCenforcement_no_column_mod_ && allEssBCs.size() > 0) {
494 enforceEssentialBC_step_2(allEssBCs);
495 }
496
497 return(0);
498}
499
500//----------------------------------------------------------------------------
501int snl_fei::LinearSystem_General::enforceEssentialBC_LinSysCore()
502{
503 fei::Matrix* matptr = matrix_.get();
504 fei::MatrixReducer* matred = dynamic_cast<fei::MatrixReducer*>(matptr);
505 if (matred != NULL) {
506 matptr = matred->getTargetMatrix().get();
507 }
508
510 dynamic_cast<fei::Matrix_Impl<LinearSystemCore>*>(matptr);
511 if (lscmatrix == 0) {
512 return(-1);
513 }
514
515 int localsize = matrixGraph_->getRowSpace()->getNumIndices_Owned();
516 fei::SharedPtr<fei::Reducer> reducer = matrixGraph_->getReducer();
517 if (matrixGraph_->getGlobalNumSlaveConstraints() > 0) {
518 localsize = reducer->getLocalReducedEqns().size();
519 }
520
521 fei::SharedPtr<fei::FillableMat> inner(new fei::FillableMat);
522 bool zeroSharedRows = false;
524 matrix.reset(new fei::Matrix_Impl<fei::FillableMat>(inner, matrixGraph_, localsize, zeroSharedRows));
525
527 matrixGraph_->getRemotelyOwnedGraphRows();
528
529 if (!BCenforcement_no_column_mod_) {
530 CHK_ERR( snl_fei::gatherRemoteEssBCs(*essBCvalues_, remoteGraph.get(), *matrix) );
531 }
532
533 unsigned numBCRows = inner->getNumRows();
534
535 if (output_stream_ != NULL && output_level_ >= fei::BRIEF_LOGS) {
536 FEI_OSTREAM& os = *output_stream_;
537 os << "#enforceEssentialBC_LinSysCore RemEssBCs to enforce: "
538 << numBCRows << FEI_ENDL;
539 }
540
541 if (numBCRows > 0 && !BCenforcement_no_column_mod_) {
542 std::vector<int*> colIndices(numBCRows);
543 std::vector<double*> coefs(numBCRows);
544 std::vector<int> colIndLengths(numBCRows);
545
546 fei::CSRMat csrmat(*inner);
547 fei::SparseRowGraph& srg = csrmat.getGraph();
548
549 int numEqns = csrmat.getNumRows();
550 int* eqns = &(srg.rowNumbers[0]);
551 int* rowOffsets = &(srg.rowOffsets[0]);
552
553 for(int i=0; i<numEqns; ++i) {
554 colIndices[i] = &(srg.packedColumnIndices[rowOffsets[i]]);
555 coefs[i] = &(csrmat.getPackedCoefs()[rowOffsets[i]]);
556 colIndLengths[i] = rowOffsets[i+1] - rowOffsets[i];
557 }
558
559 int** colInds = &colIndices[0];
560 int* colIndLens = &colIndLengths[0];
561 double** BCcoefs = &coefs[0];
562
563 if (output_stream_ != NULL && output_level_ > fei::BRIEF_LOGS) {
564 FEI_OSTREAM& os = *output_stream_;
565 for(int i=0; i<numEqns; ++i) {
566 os << "remBCeqn: " << eqns[i] << ", inds/coefs: ";
567 for(int j=0; j<colIndLens[i]; ++j) {
568 os << "("<<colInds[i][j]<<","<<BCcoefs[i][j]<<") ";
569 }
570 os << FEI_ENDL;
571 }
572 }
573
574 int errcode = lscmatrix->getMatrix()->enforceRemoteEssBCs(numEqns,
575 eqns,
576 colInds,
577 colIndLens,
578 BCcoefs);
579 if (errcode != 0) {
580 return(errcode);
581 }
582 }
583
584 int numEqns = essBCvalues_->size();
585 if (numEqns > 0) {
586 int* eqns = &(essBCvalues_->indices())[0];
587 double* bccoefs = &(essBCvalues_->coefs())[0];
588 std::vector<double> ones(numEqns, 1.0);
589
590 return(lscmatrix->getMatrix()->enforceEssentialBC(eqns, &ones[0],
591 bccoefs, numEqns));
592 }
593
594 return(0);
595}
596
597//----------------------------------------------------------------------------
598void snl_fei::LinearSystem_General::enforceEssentialBC_step_1(fei::CSVec& essBCs)
599{
600 //to enforce essential boundary conditions, we do the following:
601 //
602 // 1. for each eqn (== essBCs.indices()[n]), {
603 // put zeros in row A[eqn], but leave 1.0 on the diagonal
604 // set b[eqn] = essBCs.coefs()[n]
605 // }
606 //
607 // 2. for i in 1..numRows (i.e., all rows) {
608 // if (i in bcEqns) continue;
609 // b[i] -= A[i,eqn] * essBCs.coefs()[n]
610 // A[i,eqn] = 0.0;
611 // }
612 //
613 //It is important to note that for step 1, essBCs need only contain
614 //local eqns, but for step 2 it should contain *ALL* bc eqns.
615 //
616 //This function performs step 1.
617
618 int numEqns = essBCs.size();
619 int* eqns = &(essBCs.indices())[0];
620 double* bcCoefs = &(essBCs.coefs())[0];
621
622 std::vector<double> coefs;
623 std::vector<int> indices;
624
625 fei::SharedPtr<fei::Reducer> reducer = matrixGraph_->getReducer();
626 bool haveSlaves = reducer.get()!=NULL;
627
628 try {
629 for(int i=0; i<numEqns; i++) {
630 int eqn = eqns[i];
631
632 //if slave-constraints are present, the incoming bc-eqns are in
633 //the reduced equation space. so we actually have to translate them back
634 //to the unreduced space before passing them into the fei::Matrix object,
635 //because the fei::Matrix object has machinery to translate unreduced eqns
636 //to the reduced space.
637 //Also, our firstLocalOffset_ and lastLocalOffset_ attributes are in the
638 //unreduced space.
639 if (haveSlaves) {
640 eqn = reducer->translateFromReducedEqn(eqn);
641 }
642
643 if (eqn < firstLocalOffset_ || eqn > lastLocalOffset_) continue;
644
645 //put gamma/alpha on the rhs for this ess-BC equation.
646 double bcValue = bcCoefs[i];
647 int err = rhs_->copyIn(1, &eqn, &bcValue);
648 if (err != 0) {
649 FEI_OSTRINGSTREAM osstr;
650 osstr <<"snl_fei::LinearSystem_General::enforceEssentialBC_step_1 ERROR: "
651 << "err="<<err<<" returned from rhs_->copyIn row="<<eqn;
652 throw std::runtime_error(osstr.str());
653 }
654
655 err = getMatrixRow(matrix_.get(), eqn, coefs, indices);
656 if (err != 0 || indices.size() < 1) {
657 continue;
658 }
659
660 int rowLen = indices.size();
661 int* indPtr = &indices[0];
662
663 //first, put zeros in the row and 1.0 on the diagonal...
664 for(int j=0; j<rowLen; j++) {
665 if (indPtr[j] == eqn) coefs[j] = 1.0;
666 else coefs[j] = 0.0;
667 }
668
669 double* coefPtr = &coefs[0];
670
671 err = matrix_->copyIn(1, &eqn, rowLen, indPtr, &coefPtr);
672 if (err != 0) {
673 FEI_OSTRINGSTREAM osstr;
674 osstr <<"snl_fei::LinearSystem_General::enforceEssentialBC_step_1 ERROR: "
675 << "err="<<err<<" returned from matrix_->copyIn row="<<eqn;
676 throw std::runtime_error(osstr.str());
677 }
678 }//for i
679 }
680 catch(std::runtime_error& exc) {
681 FEI_OSTRINGSTREAM osstr;
682 osstr << "fei::LinearSystem::enforceEssentialBC: ERROR, caught exception: "
683 << exc.what();
684 throw std::runtime_error(osstr.str());
685 }
686}
687
688//----------------------------------------------------------------------------
689void snl_fei::LinearSystem_General::enforceEssentialBC_step_2(fei::CSVec& essBCs)
690{
691 //to enforce essential boundary conditions, we do the following:
692 //
693 // 1. for each eqn (== essBCs.indices()[n]), {
694 // put zeros in row A[eqn], but leave 1.0 on the diagonal
695 // set b[eqn] = essBCs.coefs()[n]
696 // }
697 //
698 // 2. for i in 1..numRows (i.e., all rows) {
699 // if (i in bcEqns) continue;
700 // b[i] -= A[i,eqn] * essBCs.coefs()[n]
701 // A[i,eqn] = 0.0;
702 // }
703 //
704 //It is important to note that for step 1, essBCs need only contain
705 //local eqns, but for step 2 it should contain *ALL* bc eqns.
706 //
707 //This function performs step 2.
708
709 int numBCeqns = essBCs.size();
710 if (numBCeqns < 1) {
711 return;
712 }
713
714 int* bcEqns = &(essBCs.indices())[0];
715 double* bcCoefs = &(essBCs.coefs())[0];
716
717 fei::SharedPtr<fei::Reducer> reducer = matrixGraph_->getReducer();
718 bool haveSlaves = reducer.get()!=NULL;
719 if (haveSlaves) {
720 for(int i=0; i<numBCeqns; ++i) {
721 bcEqns[i] = reducer->translateFromReducedEqn(bcEqns[i]);
722 }
723 }
724
725 int firstBCeqn = bcEqns[0];
726 int lastBCeqn = bcEqns[numBCeqns-1];
727
728 std::vector<double> coefs;
729 std::vector<int> indices;
730
731 int insertPoint;
732
733 int nextBCeqnOffset = 0;
734 int nextBCeqn = bcEqns[nextBCeqnOffset];
735
736 for(int i=firstLocalOffset_; i<=lastLocalOffset_; ++i) {
737 if (haveSlaves) {
738 if (reducer->isSlaveEqn(i)) continue;
739 }
740
741 bool should_continue = false;
742 if (i >= nextBCeqn) {
743 if (i == nextBCeqn) {
744 ++nextBCeqnOffset;
745 if (nextBCeqnOffset < numBCeqns) {
746 nextBCeqn = bcEqns[nextBCeqnOffset];
747 }
748 else {
749 nextBCeqn = lastLocalOffset_+1;
750 }
751
752 should_continue = true;
753 }
754 else {
755 while(nextBCeqn <= i) {
756 if (nextBCeqn == i) should_continue = true;
757 ++nextBCeqnOffset;
758 if (nextBCeqnOffset < numBCeqns) {
759 nextBCeqn = bcEqns[nextBCeqnOffset];
760 }
761 else {
762 nextBCeqn = lastLocalOffset_+1;
763 }
764 }
765 }
766 }
767
768 if (should_continue) continue;
769
770 int err = getMatrixRow(matrix_.get(), i, coefs, indices);
771 if (err != 0 || indices.size() < 1) {
772 continue;
773 }
774
775 int numIndices = indices.size();
776 int* indicesPtr = &indices[0];
777 double* coefsPtr = &coefs[0];
778 bool modifiedCoef = false;
779
780 fei::insertion_sort_with_companions(numIndices, indicesPtr, coefsPtr);
781
782 if (indicesPtr[0] > lastBCeqn || indicesPtr[numIndices-1] < firstBCeqn) {
783 continue;
784 }
785
786 double value = 0.0;
787 int offset = 0;
788
789 for(int j=0; j<numIndices; ++j) {
790 int idx = indicesPtr[j];
791 offset = fei::binarySearch(idx, bcEqns, numBCeqns,
792 insertPoint);
793 if (offset > -1) {
794 value -= bcCoefs[offset]*coefsPtr[j];
795
796 coefsPtr[j] = 0.0;
797 modifiedCoef = true;
798 }
799 }
800
801 if (modifiedCoef) {
802 err = matrix_->copyIn(1, &i, numIndices, indicesPtr, &coefsPtr);
803 if (err != 0) {
804 FEI_OSTRINGSTREAM osstr;
805 osstr <<"snl_fei::LinearSystem_General::enforceEssentialBC_step_2 ERROR: "
806 << "err="<<err<<" returned from matrix_->copyIn, row="<<i;
807 throw std::runtime_error(osstr.str());
808 }
809 }
810
811 const double fei_eps = 1.e-49;
812 if (std::abs(value) > fei_eps) {
813 rhs_->sumIn(1, &i, &value);
814
815 if (output_level_ >= fei::FULL_LOGS && output_stream_ != 0) {
816 FEI_OSTREAM& os = *output_stream_;
817 os << "enfEssBC_step2: rhs["<<i<<"] += "<<value<<FEI_ENDL;
818 }
819 }
820 }
821}
822
823//----------------------------------------------------------------------------
824int snl_fei::LinearSystem_General::getMatrixRow(fei::Matrix* matrix, int row,
825 std::vector<double>& coefs,
826 std::vector<int>& indices)
827{
828 int len = 0;
829 int err = matrix->getRowLength(row, len);
830 if (err != 0) {
831 coefs.resize(0);
832 indices.resize(0);
833 return(err);
834 }
835
836 if ((int)coefs.size() != len) {
837 coefs.resize(len);
838 }
839 if ((int)indices.size() != len) {
840 indices.resize(len);
841 }
842
843 CHK_ERR( matrix->copyOutRow(row, len, &coefs[0], &indices[0]));
844
845 return(0);
846}
847
848//----------------------------------------------------------------------------
850 const double *weights,
851 double rhsValue)
852{
853 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
854 FEI_OSTREAM& os = *output_stream_;
855 os << "loadLagrangeConstraint crID: "<<constraintID<<FEI_ENDL;
856 }
857
859 matrixGraph_->getLagrangeConstraint(constraintID);
860 if (cr == NULL) {
861 return(-1);
862 }
863
864 CHK_ERR( matrixGraph_->getConstraintConnectivityIndices(cr, iwork_) );
865
866 //Let's attach the weights to the constraint-record now.
867 std::vector<double>& cr_weights = cr->getMasterWeights();
868 cr_weights.resize(iwork_.size());
869 for(unsigned i=0; i<iwork_.size(); ++i) {
870 cr_weights.push_back(weights[i]);
871 }
872
873 fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph_->getRowSpace();
874
875 int crEqn = -1;
876 CHK_ERR( vecSpace->getGlobalIndex(cr->getIDType(),
877 cr->getConstraintID(),
878 crEqn) );
879
880 //now add the row contribution to the matrix and rhs
881 int numIndices = iwork_.size();
882 int* indicesPtr = &(iwork_[0]);
883
884 CHK_ERR( matrix_->sumIn(1, &crEqn, numIndices, indicesPtr, &weights) );
885
886 CHK_ERR( rhs_->sumIn(1, &crEqn, &rhsValue) );
887
888 //now add the column contributions to the matrix
889 for(int k=0; k<numIndices; ++k) {
890 double* thisWeight = (double*)(&(weights[k]));
891 CHK_ERR( matrix_->sumIn(1, &(indicesPtr[k]), 1, &crEqn, &thisWeight) );
892 }
893
894 return(0);
895}
896
897//----------------------------------------------------------------------------
899 const double *weights,
900 double penaltyValue,
901 double rhsValue)
902{
903 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
904 FEI_OSTREAM& os = *output_stream_;
905 os << "loadPenaltyConstraint crID: "<<constraintID<<FEI_ENDL;
906 }
907
909 matrixGraph_->getPenaltyConstraint(constraintID);
910 if (cr == NULL) {
911 return(-1);
912 }
913
914 CHK_ERR( matrixGraph_->getConstraintConnectivityIndices(cr, iwork_) );
915
916 fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph_->getRowSpace();
917
918 int numIndices = iwork_.size();
919 int* indicesPtr = &(iwork_[0]);
920
921 //now add the contributions to the matrix and rhs
922 std::vector<double> coefs(numIndices);
923 double* coefPtr = &coefs[0];
924 for(int i=0; i<numIndices; ++i) {
925 for(int j=0; j<numIndices; ++j) {
926 coefPtr[j] = weights[i]*weights[j]*penaltyValue;
927 }
928 CHK_ERR( matrix_->sumIn(1, &(indicesPtr[i]), numIndices, indicesPtr,
929 &coefPtr) );
930
931 double rhsCoef = weights[i]*penaltyValue*rhsValue;
932 CHK_ERR( rhs_->sumIn(1, &(indicesPtr[i]), &rhsCoef) );
933 }
934
935 return(0);
936}
937
virtual int loadEssentialBCs(int numIDs, const int *IDs, int idType, int fieldID, int offsetIntoField, const double *prescribedValues)
static LogManager & getLogManager()
const std::string & getOutputPath()
fei::SharedPtr< T > getMatrix()
virtual int copyOutRow(int row, int len, double *coefs, int *indices) const =0
virtual int getRowLength(int row, int &length) const =0
void reset(T *p=0)
T * get() const
std::vector< int > rowNumbers
std::vector< int > packedColumnIndices
std::vector< int > rowOffsets
virtual int copyIn(int numValues, const int *indices, const double *values, int vectorIndex=0)=0
std::vector< double > & getMasterWeights()
int loadComplete(bool applyBCs=true, bool globalAssemble=true)
LinearSystem_General(fei::SharedPtr< fei::MatrixGraph > &matrixGraph)
int loadLagrangeConstraint(int constraintID, const double *weights, double rhsValue)
bool eqnIsEssentialBC(int globalEqnIndex) const
void getConstrainedEqns(std::vector< int > &crEqns) const
int parameters(int numParams, const char *const *paramStrings)
int loadEssentialBCs(int numIDs, const int *IDs, int idType, int fieldID, int offsetIntoField, const double *prescribedValues)
int loadPenaltyConstraint(int constraintID, const double *weights, double penaltyValue, double rhsValue)
void getEssentialBCs(std::vector< int > &bcEqns, std::vector< double > &bcVals) const
void convert_ParameterSet_to_strings(const fei::ParameterSet *paramset, std::vector< std::string > &paramStrings)
Definition: fei_utils.cpp:270
fei::OutputLevel string_to_output_level(const std::string &str)
Definition: fei_utils.cpp:58
void strings_to_char_ptrs(std::vector< std::string > &stdstrings, int &numStrings, const char **&charPtrs)
Definition: fei_utils.cpp:178
int localProc(MPI_Comm comm)
int binarySearch(const T &item, const T *list, int len)
int numProcs(MPI_Comm comm)
void get_row_numbers(const FillableMat &mat, std::vector< int > &rows)
void insertion_sort_with_companions(int len, int *array, T *companions)