FEI Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
snl_fei_LinearSystem_General.cpp
Go to the documentation of this file.
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 
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 
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 {
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;
133  fei::utils::convert_ParameterSet_to_strings(&params, 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 //----------------------------------------------------------------------------
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 //----------------------------------------------------------------------------
368 void snl_fei::LinearSystem_General::getConstrainedEqns(std::vector<int>& crEqns) const
369 {
370  matrixGraph_->getConstrainedIndices(crEqns);
371 }
372 
373 //----------------------------------------------------------------------------
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 
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 //----------------------------------------------------------------------------
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 //----------------------------------------------------------------------------
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 
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 //----------------------------------------------------------------------------
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 //----------------------------------------------------------------------------
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 //----------------------------------------------------------------------------
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 
void strings_to_char_ptrs(std::vector< std::string > &stdstrings, int &numStrings, const char **&charPtrs)
Definition: fei_utils.cpp:178
virtual int getGlobalNumSlaveConstraints() const =0
virtual int loadEssentialBCs(int numIDs, const int *IDs, int idType, int fieldID, int offsetIntoField, const double *prescribedValues)
#define FEI_COUT
int loadEssentialBCs(int numIDs, const int *IDs, int idType, int fieldID, int offsetIntoField, const double *prescribedValues)
int getMatrixRow(fei::Matrix *matrix, int row, std::vector< double > &coefs, std::vector< int > &indices)
int extractDirichletBCs(fei::DirichletBCManager *bcManager, fei::SharedPtr< fei::MatrixGraph > matrixGraph, fei::CSVec *essBCvalues, bool resolveConflictRequested, bool bcs_trump_slaves)
void convert_ParameterSet_to_strings(const fei::ParameterSet *paramset, std::vector< std::string > &paramStrings)
Definition: fei_utils.cpp:270
virtual int copyIn(int numValues, const int *indices, const double *values, int vectorIndex=0)=0
const std::string & getOutputPath()
unsigned getNumRows() const
Definition: fei_CSRMat.hpp:33
void enforceEssentialBC_step_2(fei::CSVec &essBCs)
fei::OutputLevel string_to_output_level(const std::string &str)
Definition: fei_utils.cpp:58
std::vector< int > & indices()
Definition: fei_CSVec.hpp:31
virtual fei::SharedPtr< fei::Reducer > getReducer()=0
std::vector< int > rowNumbers
int parameters(int numParams, const char *const *paramStrings)
size_t size() const
Definition: fei_CSVec.hpp:36
void insertion_sort_with_companions(int len, int *array, T *companions)
int loadPenaltyConstraint(int constraintID, const double *weights, double penaltyValue, double rhsValue)
int loadComplete(bool applyBCs=true, bool globalAssemble=true)
int gatherRemoteEssBCs(fei::CSVec &essBCs, fei::SparseRowGraph *remoteGraph, fei::Matrix &matrix)
virtual int enforceEssentialBC(int *globalEqn, double *alpha, double *gamma, int len)=0
void reset(T *p=0)
void getConstrainedEqns(std::vector< int > &crEqns) const
void get_row_numbers(const FillableMat &mat, std::vector< int > &rows)
std::vector< int > packedColumnIndices
std::vector< int > rowOffsets
LinearSystem_General(fei::SharedPtr< fei::MatrixGraph > &matrixGraph)
void put_entry(CSVec &vec, int eqn, double coef)
Definition: fei_CSVec.cpp:39
int getGlobalIndex(int idType, int ID, int fieldID, int fieldOffset, int whichComponentOfField, int &globalIndex)
virtual fei::SharedPtr< fei::VectorSpace > getRowSpace()=0
void global_union(MPI_Comm comm, const fei::FillableMat &localMatrix, fei::FillableMat &globalUnionMatrix)
int binarySearch(const T &item, const T *list, int len)
std::vector< double > & getMasterWeights()
#define FEI_OSTREAM
Definition: fei_iosfwd.hpp:24
int finalizeBCEqns(fei::Matrix &matrix, bool throw_if_bc_slave_conflict=false)
int resolveConflictingCRs(fei::MatrixGraph &matrixGraph, fei::Matrix &bcEqns, const std::vector< int > &bcEqnNumbers)
void getEssentialBCs(std::vector< int > &bcEqns, std::vector< double > &bcVals) const
bool eqnIsEssentialBC(int globalEqnIndex) const
virtual int getRowLength(int row, int &length) const =0
#define ERReturn(a)
T * get() const
fei::SharedPtr< T > getMatrix()
void enforceEssentialBC_step_1(fei::CSVec &essBCs)
const char * getParam(const char *key, int numParams, const char *const *paramStrings)
fei::SharedPtr< fei::Matrix > getTargetMatrix()
int loadLagrangeConstraint(int constraintID, const double *weights, double rhsValue)
#define FEI_ENDL
SparseRowGraph & getGraph()
Definition: fei_CSRMat.hpp:27
void getGlobalIndexOffsets(std::vector< int > &globalOffsets) const
int translateFromReducedEqn(int reduced_eqn) const
void separate_BC_eqns(const fei::FillableMat &mat, std::vector< int > &bcEqns, std::vector< double > &bcVals)
std::vector< int > & getLocalReducedEqns()
int localProc(MPI_Comm comm)
static LogManager & getLogManager()
virtual int parameters(int numParams, const char *const *params)=0
virtual int enforceRemoteEssBCs(int numEqns, int *globalEqns, int **colIndices, int *colIndLen, double **coefs)=0
unsigned getNumRows() const
bool isSlaveEqn(int unreducedEqn) const
#define CHK_ERR(a)
std::vector< double > & getPackedCoefs()
Definition: fei_CSRMat.hpp:30
#define FEI_OSTRINGSTREAM
Definition: fei_sstream.hpp:32
const char * getParamValue(const char *key, int numParams, const char *const *paramStrings, char separator=' ')
virtual int copyOutRow(int row, int len, double *coefs, int *indices) const =0
int numProcs(MPI_Comm comm)
std::vector< double > & coefs()
Definition: fei_CSVec.hpp:33
int getNumIndices_Owned() const