9 #include <fei_macros.hpp> 16 #include <fei_CommUtils.hpp> 17 #include <fei_TemplateUtils.hpp> 18 #include <snl_fei_Constraint.hpp> 21 #include <fei_LibraryWrapper.hpp> 22 #include <SNL_FEI_Structure.hpp> 23 #include <fei_FiniteElementData.hpp> 24 #include <fei_Lookup.hpp> 25 #include <FEI_Implementation.hpp> 26 #include <fei_EqnCommMgr.hpp> 27 #include <fei_EqnBuffer.hpp> 28 #include <fei_NodeDatabase.hpp> 29 #include <fei_NodeCommMgr.hpp> 30 #include <fei_ProcEqns.hpp> 31 #include <fei_BlockDescriptor.hpp> 32 #include <fei_ConnectivityTable.hpp> 33 #include <snl_fei_Utils.hpp> 35 #include <fei_FEDataFilter.hpp> 38 #define fei_file "FEDataFilter.cpp" 39 #include <fei_ErrMacros.hpp> 41 #define ASSEMBLE_PUT 0 42 #define ASSEMBLE_SUM 1 49 std::vector<int>& nodeNumbers,
50 std::vector<int>& dof_ids)
52 nodeNumbers.resize(numEqns);
53 dof_ids.resize(numEqns);
55 for(
int i=0; i<numEqns; ++i) {
64 nodeNumbers[i] = nodePtr->getNodeNumber();
68 dof_ids[i] = fdmap.get_dof_id(fieldID, offset);
73 void convert_field_and_nodes_to_eqns(
const NodeDatabase& nodeDB,
74 int fieldID,
int fieldSize,
75 int numNodes,
const GlobalID* nodeIDs,
76 std::vector<int>& eqns)
78 eqns.assign(numNodes*fieldSize, -1);
81 for(
int i=0; i<numNodes; ++i) {
91 for(
int j=0; j<fieldSize; ++j) {
92 eqns[offset++] = eqn+j;
101 LibraryWrapper* wrapper,
121 masterRank_(masterRank),
122 problemStructure_(probStruct),
128 eqnCommMgr_put_(NULL),
134 constraintBlocks_(0, 16),
135 constraintNodeOffsets_(),
144 rhsIDs_.resize(numRHSs_);
147 eqnCommMgr_ = problemStructure_->getEqnCommMgr().deepCopy();
148 createEqnCommMgr_put();
150 if (wrapper_->haveFiniteElementData()) {
151 feData_ = wrapper_->getFiniteElementData();
154 fei::console_out() <<
"FEDataFilter::FEDataFilter ERROR, must be constructed with a " 155 <<
"FiniteElementData interface. Aborting." << FEI_ENDL;
157 MPI_Abort(comm_, -1);
167 char** paramStrings = NULL;
173 fei::console_out() <<
"FEDataFilter::FEDataFilter ERROR, parameters failed." << FEI_ENDL;
174 MPI_Abort(comm_, -1);
202 problemStructure_(NULL),
208 eqnCommMgr_put_(NULL),
214 constraintBlocks_(0, 16),
215 constraintNodeOffsets_(),
221 FEDataFilter::~FEDataFilter() {
228 delete eqnCommMgr_put_;
236 int FEDataFilter::initialize()
241 debugOutput(
"# initialize");
253 std::vector<int>& eqnOffsets = problemStructure_->getGlobalEqnOffsets();
254 localStartRow_ = eqnOffsets[localRank_];
255 localEndRow_ = eqnOffsets[localRank_+1]-1;
256 numGlobalEqns_ = eqnOffsets[numProcs_];
261 if (eqnCommMgr_ != NULL)
delete eqnCommMgr_;
263 if (eqnCommMgr_put_ != NULL)
delete eqnCommMgr_put_;
264 eqnCommMgr_put_ = NULL;
266 eqnCommMgr_ = problemStructure_->getEqnCommMgr().deepCopy();
267 if (eqnCommMgr_ == NULL) ERReturn(-1);
269 int err = createEqnCommMgr_put();
270 if (err != 0) ERReturn(err);
273 eqnCommMgr_->setNumRHSs(numRHSs_);
278 CHK_ERR( initLinSysCore() );
284 int FEDataFilter::createEqnCommMgr_put()
286 if (eqnCommMgr_put_ != NULL)
return(0);
288 eqnCommMgr_put_ = eqnCommMgr_->deepCopy();
289 if (eqnCommMgr_put_ == NULL) ERReturn(-1);
291 eqnCommMgr_put_->resetCoefs();
292 eqnCommMgr_put_->accumulate_ =
false;
297 int FEDataFilter::initLinSysCore()
301 int err = wrapper_->getFiniteElementData()->setLookup(*problemStructure_);
307 reducedStartRow_ = localStartRow_;
308 reducedEndRow_ = localEndRow_;
310 int numElemBlocks = problemStructure_->getNumElemBlocks();
311 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
312 NodeCommMgr& nodeCommMgr = problemStructure_->getNodeCommMgr();
315 int numRemoteNodes = nodeCommMgr.getSharedNodeIDs().size() -
316 nodeCommMgr.getLocalNodeIDs().size();
317 numNodes -= numRemoteNodes;
319 int numSharedNodes = nodeCommMgr.getNumSharedNodes();
321 std::vector<int> numElemsPerBlock(numElemBlocks);
322 std::vector<int> numNodesPerElem(numElemBlocks);
323 std::vector<int> elemMatrixSizePerBlock(numElemBlocks);
325 for(
int blk=0; blk<numElemBlocks; blk++) {
327 CHK_ERR( problemStructure_->getBlockDescriptor_index(blk, block) );
329 numElemsPerBlock[blk] = block->getNumElements();
330 numNodesPerElem[blk] = block->getNumNodesPerElement();
332 int* fieldsPerNode = block->fieldsPerNodePtr();
335 elemMatrixSizePerBlock[blk] = 0;
337 for(
int nn=0; nn<numNodesPerElem[blk]; nn++) {
338 if (fieldsPerNode[nn] <= 0) ERReturn(-1);
340 for(
int nf=0; nf<fieldsPerNode[nn]; nf++) {
341 elemMatrixSizePerBlock[blk] +=
342 problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
353 if (problemStructure_==NULL) {
354 FEI_COUT <<
"problemStructrue_ NULL"<<FEI_ENDL;
358 std::map<GlobalID,ConstraintType*>::const_iterator
359 cr_iter = problemStructure_->getPenConstRecords().begin(),
360 cr_end = problemStructure_->getPenConstRecords().end();
366 std::vector<int> numConstraintsPerBlock;
367 std::vector<int> numDofPerConstraint;
368 penCRIDs_.resize(problemStructure_->getNumPenConstRecords());
371 while(cr_iter != cr_end) {
372 penCRIDs_[counter++] = (*cr_iter).first;
376 int insertPoint = -1;
381 constraintBlocks_.insert(constraintBlocks_.begin()+insertPoint, nNodes);
382 numConstraintsPerBlock.insert(numConstraintsPerBlock.begin()+insertPoint, 1);
383 numDofPerConstraint.insert(numDofPerConstraint.begin()+insertPoint, 0);
385 if (insertPoint > 0) {
386 nodeOffset = constraintNodeOffsets_[insertPoint-1] +
387 constraintBlocks_[insertPoint-1];
389 constraintNodeOffsets_.insert(constraintNodeOffsets_.begin()+insertPoint, nodeOffset);
390 offset = insertPoint;
393 numConstraintsPerBlock[offset]++;
399 int* fieldIDs = &fieldIDsvec[0];
400 for(
int k=0; k<nNodes; ++k) {
401 int fieldSize = problemStructure_->getFieldSize(fieldIDs[k]);
402 packedFieldSizes_.insert(packedFieldSizes_.begin()+nodeOffset+k, fieldSize);
403 numDofPerConstraint[offset] += fieldSize;
409 int numBlocksTotal = numElemBlocks + constraintBlocks_.size();
410 for(
size_t i=0; i<constraintBlocks_.size(); ++i) {
411 numElemsPerBlock.push_back(numConstraintsPerBlock[i]);
412 numNodesPerElem.push_back(constraintBlocks_[i]);
413 elemMatrixSizePerBlock.push_back(numDofPerConstraint[i]);
416 int numMultCRs = problemStructure_->getNumMultConstRecords();
418 CHK_ERR( feData_->describeStructure(numBlocksTotal,
419 &numElemsPerBlock[0],
421 &elemMatrixSizePerBlock[0],
426 numRegularElems_ = 0;
427 std::vector<int> numDofPerNode;
428 std::vector<int> dof_ids;
431 for(
int i=0; i<numElemBlocks; i++) {
433 CHK_ERR( problemStructure_->getBlockDescriptor_index(i, block) );
435 if (block->getNumElements() == 0)
continue;
438 problemStructure_->getBlockConnectivity(block->getGlobalBlockID());
440 std::vector<int> cNodeList(block->getNumNodesPerElement());
442 int* fieldsPerNode = block->fieldsPerNodePtr();
445 numDofPerNode.resize(0);
446 int total_num_dof = 0;
447 for(
int nn=0; nn<numNodesPerElem[i]; nn++) {
448 if (fieldsPerNode[nn] <= 0) ERReturn(-1);
449 numDofPerNode.push_back(0);
450 int indx = numDofPerNode.size()-1;
452 for(
int nf=0; nf<fieldsPerNode[nn]; nf++) {
453 numDofPerNode[indx] += problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
455 total_num_dof += numDofPerNode[indx];
458 dof_ids.resize(total_num_dof);
460 for(
int nn=0; nn<numNodesPerElem[i]; ++nn) {
461 for(
int nf=0; nf<fieldsPerNode[nn]; ++nf) {
462 int fieldSize = problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
463 for(
int dof_offset=0; dof_offset<fieldSize; ++dof_offset) {
464 dof_ids[doffset++] = fdmap.get_dof_id(fieldIDsTable[nn][nf], dof_offset);
469 int nodesPerElement = block->getNumNodesPerElement();
472 int numElems = block->getNumElements();
473 numRegularElems_ += numElems;
474 for(
int j=0; j<numElems; j++) {
476 for(
int k=0; k<nodesPerElement; k++) {
478 cNodeList[k] = node->getNodeNumber();
481 CHK_ERR( feData_->setConnectivity(i, ctbl.elemNumbers[j],
482 block->getNumNodesPerElement(),
489 std::vector<int> nodeNumbers;
490 cr_iter = problemStructure_->getPenConstRecords().begin();
492 while(cr_iter != cr_end) {
494 std::vector<GlobalID>& nodeIDsvec = cr.
getMasters();
495 GlobalID* nodeIDs = &nodeIDsvec[0];
502 int total_num_dof = 0;
504 for(
size_t k=0; k<masterFieldIDs.size(); ++k) {
505 total_num_dof += problemStructure_->getFieldSize(masterFieldIDs[k]);
508 dof_ids.resize(total_num_dof);
510 for(
size_t k=0; k<masterFieldIDs.size(); ++k) {
511 int field_size = problemStructure_->getFieldSize(masterFieldIDs[k]);
512 for(
int dof_offset=0; dof_offset<field_size; ++dof_offset) {
513 dof_ids[doffset++] = fdmap.get_dof_id(masterFieldIDs[k], dof_offset);
517 int blockNum = numElemBlocks + index;
519 nodeNumbers.resize(nNodes);
521 for(
int k=0; k<nNodes; ++k) {
529 nodeNumbers[k] = node->getNodeNumber();
533 int offset = constraintNodeOffsets_[index];
534 CHK_ERR( feData_->setConnectivity(blockNum, numRegularElems_+i++, nNodes, &nodeNumbers[0], &packedFieldSizes_[offset], &dof_ids[0]) );
539 catch(std::runtime_error& exc) {
548 int FEDataFilter::resetSystem(
double s)
553 if (Filter::logStream() != NULL) {
554 (*logStream()) <<
"FEI: resetSystem" << FEI_ENDL << s << FEI_ENDL;
557 CHK_ERR( feData_->reset() );
559 debugOutput(
"#FEDataFilter leaving resetSystem");
565 int FEDataFilter::deleteMultCRs()
568 debugOutput(
"#FEDataFilter::deleteMultCRs");
570 int err = feData_->deleteConstraints();
572 debugOutput(
"#FEDataFilter leaving deleteMultCRs");
578 int FEDataFilter::resetTheMatrix(
double s)
586 int FEDataFilter::resetTheRHSVector(
double s)
594 int FEDataFilter::resetMatrix(
double s)
600 debugOutput(
"FEI: resetMatrix");
602 CHK_ERR( resetTheMatrix(s) );
604 eqnCommMgr_->resetCoefs();
606 debugOutput(
"#FEDataFilter leaving resetMatrix");
612 int FEDataFilter::resetRHSVector(
double s)
618 debugOutput(
"FEI: resetRHSVector");
620 CHK_ERR( resetTheRHSVector(s) );
622 eqnCommMgr_->resetCoefs();
624 debugOutput(
"#FEDataFilter leaving resetRHSVector");
630 int FEDataFilter::resetInitialGuess(
double s)
635 if (Filter::logStream() != NULL) {
636 FEI_OSTREAM& os = *logStream();
637 os <<
"FEI: resetInitialGuess" << FEI_ENDL;
638 os <<
"#value to which initial guess is to be set" << FEI_ENDL;
645 debugOutput(
"#FEDataFilter leaving resetInitialGuess");
651 int FEDataFilter::loadNodeBCs(
int numNodes,
652 const GlobalID *nodeIDs,
654 const int* offsetsIntoField,
655 const double* prescribedValues)
660 int size = problemStructure_->getFieldSize(fieldID);
662 fei::console_out() <<
"FEI Warning: loadNodeBCs called for fieldID "<<fieldID
663 <<
", which was defined with size "<<size<<
" (should be positive)."<<FEI_ENDL;
667 if (Filter::logStream() != NULL) {
668 (*logStream())<<
"FEI: loadNodeBCs"<<FEI_ENDL
669 <<
"#num-nodes"<<FEI_ENDL<<numNodes<<FEI_ENDL
670 <<
"#fieldID"<<FEI_ENDL<<fieldID<<FEI_ENDL
671 <<
"#field-size"<<FEI_ENDL<<size<<FEI_ENDL;
672 (*logStream())<<
"#following lines: nodeID offsetIntoField value "<<FEI_ENDL;
674 for(
int j=0; j<numNodes; j++) {
675 int nodeID = nodeIDs[j];
676 (*logStream())<<nodeID<<
" ";
677 (*logStream())<< offsetsIntoField[j]<<
" ";
678 (*logStream())<< prescribedValues[j]<<FEI_ENDL;
682 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
684 std::vector<int> essEqns(numNodes);
685 std::vector<double> alpha(numNodes);
686 std::vector<double> gamma(numNodes);
688 for(
int i=0; i<numNodes; ++i) {
692 fei::console_out() <<
"fei_FEDataFilter::loadNodeBCs ERROR, node " << nodeIDs[i]
693 <<
" not found." << FEI_ENDL;
702 essEqns[i] = eqn + offsetsIntoField[i];
703 gamma[i] = prescribedValues[i];
707 if (essEqns.size() > 0) {
708 CHK_ERR( enforceEssentialBCs(&essEqns[0], &alpha[0],
709 &gamma[0], essEqns.size()) );
716 int FEDataFilter::loadElemBCs(
int numElems,
717 const GlobalID *elemIDs,
719 const double *
const *alpha,
720 const double *
const *beta,
721 const double *
const *gamma)
727 void FEDataFilter::allocElemStuff()
729 int nb = problemStructure_->getNumElemBlocks();
731 for(
int i=0; i<nb; i++) {
733 int err = problemStructure_->getBlockDescriptor_index(i, block);
734 if (err) voidERReturn;
736 int numEqns = block->getNumEqnsPerElement();
737 if (maxElemRows_ < numEqns) maxElemRows_ = numEqns;
740 eStiff_ =
new double*[maxElemRows_];
741 eStiff1D_ =
new double[maxElemRows_*maxElemRows_];
743 if (eStiff_ == NULL || eStiff1D_ == NULL) voidERReturn
745 for(
int r=0; r<maxElemRows_; r++) {
746 eStiff_[r] = eStiff1D_ + r*maxElemRows_;
749 eLoad_ =
new double[maxElemRows_];
751 if (eLoad_ == NULL) voidERReturn
755 int FEDataFilter::sumInElem(GlobalID elemBlockID,
757 const GlobalID* elemConn,
758 const double*
const* elemStiffness,
759 const double* elemLoad,
762 if (Filter::logStream() != NULL) {
763 (*logStream()) <<
"FEI: sumInElem" << FEI_ENDL <<
"# elemBlockID " << FEI_ENDL
764 << static_cast<int>(elemBlockID) << FEI_ENDL
765 <<
"# elemID " << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
767 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
768 int numNodes = block->getNumNodesPerElement();
769 (*logStream()) <<
"#num-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
770 (*logStream()) <<
"#connected nodes" << FEI_ENDL;
771 for(
int i=0; i<numNodes; ++i) {
772 GlobalID nodeID = elemConn[i];
773 (*logStream())<<static_cast<int>(nodeID)<<
" ";
775 (*logStream())<<FEI_ENDL;
778 return(generalElemInput(elemBlockID, elemID, elemConn, elemStiffness,
779 elemLoad, elemFormat));
783 int FEDataFilter::sumInElemMatrix(GlobalID elemBlockID,
785 const GlobalID* elemConn,
786 const double*
const* elemStiffness,
789 if (Filter::logStream() != NULL) {
790 (*logStream()) <<
"FEI: sumInElemMatrix"<<FEI_ENDL
791 <<
"#elemBlockID" << FEI_ENDL << static_cast<int>(elemBlockID)
792 <<
"# elemID" << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
794 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
795 int numNodes = block->getNumNodesPerElement();
796 (*logStream()) <<
"#num-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
797 (*logStream()) <<
"#connected nodes" << FEI_ENDL;
798 for(
int i=0; i<numNodes; ++i) {
799 GlobalID nodeID = elemConn[i];
800 (*logStream())<<static_cast<int>(nodeID)<<
" ";
802 (*logStream())<<FEI_ENDL;
805 return(generalElemInput(elemBlockID, elemID, elemConn, elemStiffness,
810 int FEDataFilter::sumInElemRHS(GlobalID elemBlockID,
812 const GlobalID* elemConn,
813 const double* elemLoad)
815 if (Filter::logStream() != NULL) {
816 (*logStream()) <<
"FEI: sumInElemRHS"<<FEI_ENDL<<
"# elemBlockID " << FEI_ENDL
817 <<static_cast<int>(elemBlockID)
818 <<
"# elemID " << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
820 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
821 int numNodes = block->getNumNodesPerElement();
822 (*logStream()) <<
"#num-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
823 (*logStream()) <<
"#connected nodes" << FEI_ENDL;
824 for(
int i=0; i<numNodes; ++i) {
825 GlobalID nodeID = elemConn[i];
826 (*logStream())<<static_cast<int>(nodeID)<<
" ";
828 (*logStream())<<FEI_ENDL;
831 return(generalElemInput(elemBlockID, elemID, elemConn, NULL,
836 int FEDataFilter::generalElemInput(GlobalID elemBlockID,
838 const GlobalID* elemConn,
839 const double*
const* elemStiffness,
840 const double* elemLoad,
844 return(generalElemInput(elemBlockID, elemID, elemStiffness, elemLoad,
849 int FEDataFilter::generalElemInput(GlobalID elemBlockID,
851 const double*
const* elemStiffness,
852 const double* elemLoad,
858 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
862 if (maxElemRows_ <= 0) allocElemStuff();
864 int numElemRows = block->getNumEqnsPerElement();
869 const double*
const* stiff = NULL;
870 if (elemStiffness != NULL) stiff = elemStiffness;
872 const double* load = NULL;
873 if (elemLoad != NULL) load = elemLoad;
878 if (elemFormat != FEI_DENSE_ROW && stiff != NULL) {
879 Filter::copyStiffness(stiff, numElemRows, elemFormat, eStiff_);
883 if (stiff != NULL || load != NULL) newData_ =
true;
885 if (Filter::logStream() != NULL) {
888 <<
"#numElemRows"<< FEI_ENDL << numElemRows << FEI_ENDL
889 <<
"#elem-stiff (after being copied into dense-row format)" 891 for(
int i=0; i<numElemRows; i++) {
892 const double* stiff_i = stiff[i];
893 for(
int j=0; j<numElemRows; j++) {
894 (*logStream()) << stiff_i[j] <<
" ";
896 (*logStream()) << FEI_ENDL;
901 (*logStream()) <<
"#elem-load" << FEI_ENDL;
902 for(
int i=0; i<numElemRows; i++) {
903 (*logStream()) << load[i] <<
" ";
905 (*logStream())<<FEI_ENDL;
912 int blockNumber = problemStructure_->getIndexOfBlock(elemBlockID);
915 getBlockConnectivity(elemBlockID);
917 std::map<GlobalID,int>::iterator
918 iter = connTable.elemIDs.find(elemID);
919 if (iter == connTable.elemIDs.end()) {
925 int elemIndex = iter->second;
927 int elemNumber = connTable.elemNumbers[elemIndex];
929 int numNodes = block->getNumNodesPerElement();
930 int* fieldsPerNode = block->fieldsPerNodePtr();
933 int numDistinctFields = block->getNumDistinctFields();
936 int total_num_dofs = 0;
937 if (numDistinctFields == 1) {
938 fieldSize = problemStructure_->getFieldSize(fieldIDsTable[0][0]);
939 for(
int i=0; i<numNodes; ++i) {
940 total_num_dofs += fieldSize*fieldsPerNode[i];
942 dof_id = fdmap.get_dof_id(fieldIDsTable[0][0], 0);
945 for(
int i=0; i<numNodes; ++i) {
946 for(
int nf=0; nf<fieldsPerNode[i]; ++nf) {
947 total_num_dofs += problemStructure_->getFieldSize(fieldIDsTable[i][nf]);
952 static std::vector<int> iwork;
953 iwork.resize(2*numNodes+total_num_dofs);
955 int* dofsPerNode = &iwork[0];
956 int* nodeNumbers = dofsPerNode+numNodes;
957 int* dof_ids = nodeNumbers+numNodes;
959 for(
int i=0; i<numNodes; ++i) {
965 &((*connTable.elem_conn_ptrs)[elemIndex*numNodes]);
968 for(
int nn=0; nn<numNodes; nn++) {
970 nodeNumbers[nn] = node->getNodeNumber();
972 if (numDistinctFields == 1) {
973 for(
int nf=0; nf<fieldsPerNode[nn]; nf++) {
974 dofsPerNode[nn] += fieldSize;
975 for(
int dof_offset=0; dof_offset<fieldSize; ++dof_offset) {
976 dof_ids[doffset++] = dof_id;
981 for(
int nf=0; nf<fieldsPerNode[nn]; nf++) {
982 int fieldSize = problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
983 int dof_id = fdmap.get_dof_id(fieldIDsTable[nn][nf], 0);
984 dofsPerNode[nn] += fieldSize;
985 for(
int dof_offset=0; dof_offset<fieldSize; ++dof_offset) {
986 dof_ids[doffset++] = dof_id + dof_offset;
993 CHK_ERR( feData_->setElemMatrix(blockNumber, elemNumber, numNodes,
994 nodeNumbers, dofsPerNode, dof_ids, stiff) );
998 CHK_ERR( feData_->setElemVector(blockNumber, elemNumber, numNodes,
999 nodeNumbers, dofsPerNode, dof_ids, load) );
1002 return(FEI_SUCCESS);
1006 int FEDataFilter::putIntoRHS(
int IDType,
1009 const GlobalID* IDs,
1010 const double* rhsEntries)
1012 int fieldSize = problemStructure_->getFieldSize(fieldID);
1014 rowIndices_.resize(fieldSize*numIDs);
1017 CHK_ERR( problemStructure_->getEqnNumbers(numIDs, IDs, IDType, fieldID,
1020 if (checkNumEqns != numIDs*fieldSize) {
1024 CHK_ERR( exchangeRemoteEquations() );
1026 CHK_ERR(assembleRHS(rowIndices_.size(), &rowIndices_[0], rhsEntries, ASSEMBLE_PUT));
1032 int FEDataFilter::sumIntoRHS(
int IDType,
1035 const GlobalID* IDs,
1036 const double* rhsEntries)
1038 int fieldSize = problemStructure_->getFieldSize(fieldID);
1040 rowIndices_.resize(fieldSize*numIDs);
1043 CHK_ERR( problemStructure_->getEqnNumbers(numIDs, IDs, IDType, fieldID,
1046 if (checkNumEqns != numIDs*fieldSize) {
1050 CHK_ERR(assembleRHS(rowIndices_.size(), &rowIndices_[0], rhsEntries, ASSEMBLE_SUM));
1056 int FEDataFilter::sumIntoMatrixDiagonal(
int IDType,
1059 const GlobalID* IDs,
1060 const double* coefficients)
1062 int fieldSize = problemStructure_->getFieldSize(fieldID);
1063 const NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1065 std::vector<int> eqns;
1066 convert_field_and_nodes_to_eqns(nodeDB, fieldID, fieldSize, numIDs, IDs, eqns);
1068 std::vector<int> nodeNumbers, dof_ids;
1069 convert_eqns_to_nodenumbers_and_dof_ids(problemStructure_->getFieldDofMap(),
1070 nodeDB, eqns.size(), &eqns[0],
1071 nodeNumbers, dof_ids);
1073 std::vector<int> ones(nodeNumbers.size(), 1);
1075 CHK_ERR( feData_->sumIntoMatrix(nodeNumbers.size(), &nodeNumbers[0], &dof_ids[0],
1076 &ones[0], &nodeNumbers[0], &dof_ids[0], coefficients));
1081 int FEDataFilter::enforceEssentialBCs(
const int* eqns,
1082 const double* alpha,
1083 const double* gamma,
1086 std::vector<double> values;
1087 std::vector<int> nodeNumbers;
1088 std::vector<int> dof_ids;
1091 for(
int i=0; i<numEqns; i++) {
1092 int reducedEqn = -1;
1093 bool isSlave = problemStructure_->
1094 translateToReducedEqn(eqns[i], reducedEqn);
1095 if (isSlave)
continue;
1097 int nodeNumber = problemStructure_->getAssociatedNodeNumber(eqns[i]);
1099 nodeNumbers.push_back(nodeNumber);
1102 CHK_ERR( problemStructure_->getNodeDatabase().
1103 getNodeWithNumber(nodeNumber, node));
1105 int fieldID, offset;
1107 dof_ids.push_back( fdmap.get_dof_id(fieldID, offset) );
1109 values.push_back(gamma[i]/alpha[i]);
1112 CHK_ERR( feData_->setDirichletBCs(nodeNumbers.size(),
1113 &nodeNumbers[0], &dof_ids[0], &values[0]) );
1117 return(FEI_SUCCESS);
1121 int FEDataFilter::loadFEDataMultCR(
int CRID,
1123 const GlobalID* CRNodes,
1124 const int* CRFields,
1125 const double* CRWeights,
1128 if (Filter::logStream() != NULL) {
1129 FEI_OSTREAM& os = *logStream();
1130 os<<
"FEI: loadCRMult"<<FEI_ENDL;
1131 os<<
"#num-nodes"<<FEI_ENDL<<numCRNodes<<FEI_ENDL;
1132 os<<
"#CRNodes:"<<FEI_ENDL;
1134 for(i=0; i<numCRNodes; ++i) {
1135 GlobalID nodeID = CRNodes[i];
1136 os << static_cast<int>(nodeID) <<
" ";
1138 os << FEI_ENDL <<
"#fields:"<<FEI_ENDL;
1139 for(i=0; i<numCRNodes; ++i) os << CRFields[i] <<
" ";
1140 os << FEI_ENDL <<
"#field-sizes:"<<FEI_ENDL;
1141 for(i=0; i<numCRNodes; ++i) {
1142 int size = problemStructure_->getFieldSize(CRFields[i]);
1145 os << FEI_ENDL<<
"#weights:"<<FEI_ENDL;
1147 for(i=0; i<numCRNodes; ++i) {
1148 int size = problemStructure_->getFieldSize(CRFields[i]);
1149 for(
int j=0; j<size; ++j) {
1150 os << CRWeights[offset++] <<
" ";
1153 os << FEI_ENDL<<
"#CRValue:"<<FEI_ENDL<<CRValue<<FEI_ENDL;
1156 if (numCRNodes <= 0)
return(0);
1158 std::vector<int> nodeNumbers;
1159 std::vector<int> dof_ids;
1160 std::vector<double> weights;
1162 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1165 double fei_min = std::numeric_limits<double>::min();
1168 for(
int i=0; i<numCRNodes; i++) {
1174 if (!hasField) ERReturn(-1);
1176 int fieldSize = problemStructure_->getFieldSize(CRFields[i]);
1177 int dof_id = fdmap.get_dof_id(CRFields[i], 0);
1179 for(
int f=0; f<fieldSize; f++) {
1180 double weight = CRWeights[offset++];
1181 if (std::abs(weight) > fei_min) {
1182 nodeNumbers.push_back(node->getNodeNumber());
1183 dof_ids.push_back(dof_id+f);
1184 weights.push_back(weight);
1189 CHK_ERR( feData_->setMultiplierCR(CRID, nodeNumbers.size(),
1190 &nodeNumbers[0], &dof_ids[0],
1191 &weights[0], CRValue) );
1198 int FEDataFilter::loadFEDataPenCR(
int CRID,
1200 const GlobalID* CRNodes,
1201 const int* CRFields,
1202 const double* CRWeights,
1206 if (numCRNodes <= 0)
return(0);
1208 std::vector<int> nodeNumbers;
1209 std::vector<int> dofsPerNode;
1210 std::vector<int> dof_ids;
1211 std::vector<double> weights;
1213 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1217 for(
int i=0; i<numCRNodes; i++) {
1220 if(node == NULL)
continue;
1225 if (!hasField)
continue;
1227 int fieldSize = problemStructure_->getFieldSize(CRFields[i]);
1229 nodeNumbers.push_back(node->getNodeNumber());
1230 dofsPerNode.push_back(fieldSize);
1232 for(
int f=0; f<fieldSize; f++) {
1233 dof_ids.push_back(fdmap.get_dof_id(CRFields[i], f));
1234 double weight = CRWeights[offset++];
1235 weights.push_back(weight);
1239 std::vector<double*> matrixCoefs(weights.size());
1240 std::vector<double> rhsCoefs(weights.size());
1242 for(
size_t i=0; i<weights.size(); ++i) {
1243 double* coefPtr =
new double[weights.size()];
1244 for(
size_t j=0; j<weights.size(); ++j) {
1245 coefPtr[j] = weights[i]*weights[j]*penValue;
1247 matrixCoefs[i] = coefPtr;
1248 rhsCoefs[i] = weights[i]*penValue*CRValue;
1255 int blockNum = problemStructure_->getNumElemBlocks() + index;
1256 int elemNum = numRegularElems_ + crIndex;
1258 CHK_ERR( feData_->setElemMatrix(blockNum, elemNum,
1265 CHK_ERR( feData_->setElemVector(blockNum, elemNum, nodeNumbers.size(),
1266 &nodeNumbers[0], &dofsPerNode[0], &dof_ids[0], &rhsCoefs[0]) );
1270 for(
size_t i=0; i<weights.size(); ++i) {
1271 delete [] matrixCoefs[i];
1278 int FEDataFilter::loadCRMult(
int CRID,
1280 const GlobalID* CRNodes,
1281 const int* CRFields,
1282 const double* CRWeights,
1290 CHK_ERR( loadFEDataMultCR(CRID, numCRNodes, CRNodes, CRFields, CRWeights,
1293 return(FEI_SUCCESS);
1297 int FEDataFilter::loadCRPen(
int CRID,
1299 const GlobalID* CRNodes,
1300 const int* CRFields,
1301 const double* CRWeights,
1309 debugOutput(
"FEI: loadCRPen");
1312 CHK_ERR( loadFEDataPenCR(CRID, numCRNodes, CRNodes, CRFields, CRWeights,
1313 CRValue, penValue) );
1315 return(FEI_SUCCESS);
1319 int FEDataFilter::parameters(
int numParams,
const char *
const* paramStrings)
1325 if (numParams == 0 || paramStrings == NULL) {
1326 debugOutput(
"#FEDataFilter::parameters --- no parameters.");
1330 snl_fei::getIntParamValue(
"outputLevel",numParams, paramStrings,
1333 snl_fei::getIntParamValue(
"internalFei",numParams,paramStrings,
1336 if (Filter::logStream() != NULL) {
1337 (*logStream())<<
"#FEDataFilter::parameters"<<FEI_ENDL
1338 <<
"# --- numParams: "<< numParams<<FEI_ENDL;
1339 for(
int i=0; i<numParams; i++){
1340 (*logStream())<<
"#------ paramStrings["<<i<<
"]: " 1341 <<paramStrings[i]<<FEI_ENDL;
1346 CHK_ERR( Filter::parameters(numParams, paramStrings) );
1348 debugOutput(
"#FEDataFilter leaving parameters function");
1350 return(FEI_SUCCESS);
1354 int FEDataFilter::loadComplete()
1356 debugOutput(
"#FEDataFilter calling FEData matrixLoadComplete");
1358 CHK_ERR( feData_->loadComplete() );
1366 int FEDataFilter::residualNorm(
int whichNorm,
int numFields,
1367 int* fieldIDs,
double* norms,
double& residTime)
1373 debugOutput(
"FEI: residualNorm");
1375 CHK_ERR( loadComplete() );
1379 int fdbNumFields = problemStructure_->getNumFields();
1380 const int* fdbFieldIDs = problemStructure_->getFieldIDsPtr();
1388 while(offset < numFields && i < fdbNumFields) {
1389 if (fdbFieldIDs[i] >= 0) {
1390 fieldIDs[offset++] = fdbFieldIDs[i];
1394 for(i=0; i<numFields; ++i) {
1400 for(i=offset; i<numFields; ++i) {
1404 return(FEI_SUCCESS);
1408 int FEDataFilter::formResidual(
double* residValues,
int numLocalEqns)
1411 return(FEI_SUCCESS);
1415 int FEDataFilter::solve(
int& status,
double& sTime) {
1417 debugOutput(
"FEI: solve");
1419 CHK_ERR( loadComplete() );
1421 debugOutput(
"#FEDataFilter in solve, calling launchSolver...");
1423 double start = MPI_Wtime();
1425 CHK_ERR( feData_->launchSolver(status, iterations_) );
1427 sTime = MPI_Wtime() - start;
1429 debugOutput(
"#FEDataFilter... back from solver");
1433 CHK_ERR( unpackSolution() );
1435 debugOutput(
"#FEDataFilter leaving solve");
1437 if (status != 0)
return(1);
1438 else return(FEI_SUCCESS);
1442 int FEDataFilter::setNumRHSVectors(
int numRHSs,
int* rhsIDs){
1445 fei::console_out() <<
"FEDataFilter::setNumRHSVectors: ERROR, numRHSs < 0." << FEI_ENDL;
1451 rhsIDs_.resize(numRHSs_);
1452 for(
int i=0; i<numRHSs_; i++) rhsIDs_[i] = rhsIDs[i];
1455 eqnCommMgr_->setNumRHSs(numRHSs_);
1457 return(FEI_SUCCESS);
1461 int FEDataFilter::setCurrentRHS(
int rhsID)
1463 std::vector<int>::iterator iter =
1464 std::find( rhsIDs_.begin(), rhsIDs_.end(), rhsID);
1466 if (iter == rhsIDs_.end()) ERReturn(-1)
1468 currentRHS_ = iter - rhsIDs_.begin();
1470 return(FEI_SUCCESS);
1474 int FEDataFilter::giveToMatrix(
int numPtRows,
const int* ptRows,
1475 int numPtCols,
const int* ptCols,
1476 const double*
const* values,
int mode)
1481 std::vector<int> rowNodeNumbers, row_dof_ids, colNodeNumbers, col_dof_ids;
1482 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1490 for(i=0; i<numPtRows; i++) {
1491 int nodeNumber = problemStructure_->getAssociatedNodeNumber(ptRows[i]);
1492 if (nodeNumber < 0) ERReturn(-1);
1495 int fieldID, offset;
1496 node->
getFieldID(ptRows[i], fieldID, offset);
1498 rowNodeNumbers.push_back(nodeNumber);
1499 row_dof_ids.push_back(fdmap.get_dof_id(fieldID, offset));
1502 for(i=0; i<numPtCols; i++) {
1503 int nodeNumber = problemStructure_->getAssociatedNodeNumber(ptCols[i]);
1504 if (nodeNumber < 0) ERReturn(-1);
1507 int fieldID, offset;
1508 node->
getFieldID(ptCols[i], fieldID, offset);
1510 colNodeNumbers.push_back(nodeNumber);
1511 col_dof_ids.push_back(fdmap.get_dof_id(fieldID, offset));
1518 int len = numPtRows*numPtCols;
1519 std::vector<int> allColNodeNumbers(len), all_col_dof_ids(len);
1520 std::vector<double> allValues(len);
1523 for(i=0; i<numPtRows; i++) {
1524 for(
int j=0; j<numPtCols; j++) {
1525 allColNodeNumbers[offset] = colNodeNumbers[j];
1526 all_col_dof_ids[offset] = col_dof_ids[j];
1527 allValues[offset++] = values[i][j];
1533 std::vector<int> numColsPerRow(numPtRows, numPtCols);
1538 CHK_ERR( feData_->sumIntoMatrix(numPtRows,
1542 &allColNodeNumbers[0],
1543 &all_col_dof_ids[0],
1546 return(FEI_SUCCESS);
1550 int FEDataFilter::giveToLocalReducedMatrix(
int numPtRows,
const int* ptRows,
1551 int numPtCols,
const int* ptCols,
1552 const double*
const* values,
int mode)
1557 std::vector<int> rowNodeNumbers, row_dof_ids, colNodeNumbers, col_dof_ids;
1558 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1566 for(i=0; i<numPtRows; i++) {
1567 int nodeNumber = problemStructure_->getAssociatedNodeNumber(ptRows[i]);
1568 if (nodeNumber < 0) ERReturn(-1);
1571 int fieldID, offset;
1572 node->
getFieldID(ptRows[i], fieldID, offset);
1574 rowNodeNumbers.push_back(nodeNumber);
1575 row_dof_ids.push_back(fdmap.get_dof_id(fieldID, offset));
1578 for(i=0; i<numPtCols; i++) {
1579 int nodeNumber = problemStructure_->getAssociatedNodeNumber(ptCols[i]);
1580 if (nodeNumber < 0) ERReturn(-1);
1583 int fieldID, offset;
1584 node->
getFieldID(ptCols[i], fieldID, offset);
1586 colNodeNumbers.push_back(nodeNumber);
1587 col_dof_ids.push_back(fdmap.get_dof_id(fieldID, offset));
1594 int len = numPtRows*numPtCols;
1595 std::vector<int> allColNodeNumbers(len), all_col_dof_ids(len);
1596 std::vector<double> allValues(len);
1599 for(i=0; i<numPtRows; i++) {
1600 for(
int j=0; j<numPtCols; j++) {
1601 allColNodeNumbers[offset] = colNodeNumbers[j];
1602 all_col_dof_ids[offset] = col_dof_ids[j];
1603 allValues[offset++] = values[i][j];
1609 std::vector<int> numColsPerRow(numPtRows, numPtCols);
1614 CHK_ERR( feData_->sumIntoMatrix(numPtRows,
1618 &allColNodeNumbers[0],
1619 &all_col_dof_ids[0],
1622 return(FEI_SUCCESS);
1626 int FEDataFilter::getFromMatrix(
int numPtRows,
const int* ptRows,
1627 const int* rowColOffsets,
const int* ptCols,
1628 int numColsPerRow,
double** values)
1647 int FEDataFilter::giveToRHS(
int num,
const double* values,
1648 const int* indices,
int mode)
1650 std::vector<int> workspace(num*2);
1651 int* rowNodeNumbers = &workspace[0];
1652 int* row_dof_ids = rowNodeNumbers+num;
1653 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1656 for(
int i=0; i<num; ++i) {
1660 rowNodeNumbers[i] = -1;
1661 row_dof_ids[i] = -1;
1665 rowNodeNumbers[i] = nodeptr->getNodeNumber();
1667 int fieldID, offset;
1668 nodeptr->
getFieldID(indices[i], fieldID, offset);
1670 row_dof_ids[i] = fdmap.get_dof_id(fieldID, offset);
1673 if (mode == ASSEMBLE_SUM) {
1674 CHK_ERR( feData_->sumIntoRHSVector(num,
1680 CHK_ERR( feData_->putIntoRHSVector(num,
1686 return(FEI_SUCCESS);
1690 int FEDataFilter::giveToLocalReducedRHS(
int num,
const double* values,
1691 const int* indices,
int mode)
1693 std::vector<int> rowNodeNumbers, row_dof_ids;
1694 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1697 for(
int i=0; i<num; i++) {
1698 int nodeNumber = problemStructure_->getAssociatedNodeNumber(indices[i]);
1699 if (nodeNumber < 0) ERReturn(-1);
1704 int fieldID, offset;
1705 node->
getFieldID(indices[i], fieldID, offset);
1707 rowNodeNumbers.push_back(nodeNumber);
1708 row_dof_ids.push_back(fdmap.get_dof_id(fieldID, offset));
1711 if (mode == ASSEMBLE_SUM) {
1712 CHK_ERR( feData_->sumIntoRHSVector(rowNodeNumbers.size(),
1714 &row_dof_ids[0], values) );
1717 CHK_ERR( feData_->putIntoRHSVector(rowNodeNumbers.size(),
1719 &row_dof_ids[0], values) );
1722 return(FEI_SUCCESS);
1726 int FEDataFilter::getFromRHS(
int num,
double* values,
const int* indices)
1732 int FEDataFilter::getEqnSolnEntry(
int eqnNumber,
double& solnValue)
1741 if (localStartRow_ > eqnNumber || eqnNumber > localEndRow_) {
1743 CHK_ERR( getSharedRemoteSolnEntry(eqnNumber, solnValue) );
1747 CHK_ERR( getReducedSolnEntry( eqnNumber, solnValue ) );
1753 int FEDataFilter::getSharedRemoteSolnEntry(
int eqnNumber,
double& solnValue)
1755 std::vector<int>& remoteEqnNumbers = eqnCommMgr_->sendEqnNumbersPtr();
1756 double* remoteSoln = eqnCommMgr_->sendEqnSolnPtr();
1760 fei::console_out() <<
"FEDataFilter::getSharedRemoteSolnEntry: ERROR, eqn " 1761 << eqnNumber <<
" not found." << FEI_ENDL;
1764 solnValue = remoteSoln[index];
1769 int FEDataFilter::getReducedSolnEntry(
int eqnNumber,
double& solnValue)
1776 int nodeNumber = problemStructure_->getAssociatedNodeNumber(eqnNumber);
1781 if (nodeNumber < 0) {solnValue = -999.99;
return(FEI_SUCCESS);}
1784 problemStructure_->getNodeDatabase().getNodeWithNumber(nodeNumber, node);
1794 int eqn = problemStructure_->translateFromReducedEqn(eqnNumber);
1795 int fieldID, offset;
1797 int dof_id = problemStructure_->getFieldDofMap().get_dof_id(fieldID, offset);
1799 bool fetiHasNode =
true;
1800 GlobalID nodeID = node->getGlobalNodeID();
1801 NodeCommMgr& nodeCommMgr = problemStructure_->getNodeCommMgr();
1802 std::vector<GlobalID>& shNodeIDs = nodeCommMgr.getSharedNodeIDs();
1805 if (!(problemStructure_->isInLocalElement(nodeNumber)) ) fetiHasNode =
false;
1809 int err = feData_->getSolnEntry(nodeNumber, dof_id, solnValue);
1811 fei::console_out() <<
"FEDataFilter::getReducedSolnEntry: nodeNumber " << nodeNumber
1812 <<
" (nodeID " << node->getGlobalNodeID() <<
"), dof_id "<<dof_id
1813 <<
" couldn't be obtained from FETI on proc " << localRank_ << FEI_ENDL;
1818 return(FEI_SUCCESS);
1822 int FEDataFilter::unpackSolution()
1830 if (Filter::logStream() != NULL) {
1831 (*logStream())<<
"# entering unpackSolution, outputLevel: " 1832 <<outputLevel_<<FEI_ENDL;
1841 int numRecvEqns = eqnCommMgr_->getNumLocalEqns();
1842 std::vector<int>& recvEqnNumbers = eqnCommMgr_->localEqnNumbers();
1844 for(
int i=0; i<numRecvEqns; i++) {
1845 int eqn = recvEqnNumbers[i];
1847 if ((reducedStartRow_ > eqn) || (reducedEndRow_ < eqn)) {
1848 fei::console_out() <<
"FEDataFilter::unpackSolution: ERROR, 'recv' eqn (" << eqn
1849 <<
") out of local range." << FEI_ENDL;
1850 MPI_Abort(comm_, -1);
1853 double solnValue = 0.0;
1855 CHK_ERR( getReducedSolnEntry(eqn, solnValue) );
1857 eqnCommMgr_->addSolnValues(&eqn, &solnValue, 1);
1860 eqnCommMgr_->exchangeSoln();
1862 debugOutput(
"#FEDataFilter leaving unpackSolution");
1863 return(FEI_SUCCESS);
1867 void FEDataFilter:: setEqnCommMgr(
EqnCommMgr* eqnCommMgr)
1870 eqnCommMgr_ = eqnCommMgr;
1874 int FEDataFilter::getBlockNodeSolution(GlobalID elemBlockID,
1876 const GlobalID *nodeIDs,
1880 debugOutput(
"FEI: getBlockNodeSolution");
1882 int numActiveNodes = problemStructure_->getNumActiveNodes();
1883 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1885 if (numActiveNodes <= 0)
return(0);
1887 int numSolnParams = 0;
1890 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
1896 for(
int i=0; i<numActiveNodes; i++) {
1900 if (offset == numNodes)
break;
1902 GlobalID nodeID = nodeIDs[offset];
1905 offsets[offset++] = numSolnParams;
1912 if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
1926 int numFields = node->getNumFields();
1927 const int* fieldIDs = node->getFieldIDList();
1929 for(
int j=0; j<numFields; j++) {
1930 if (block->containsField(fieldIDs[j])) {
1931 int size = problemStructure_->getFieldSize(fieldIDs[j]);
1940 for(
int k=0; k<size; k++) {
1941 CHK_ERR( getEqnSolnEntry(thisEqn+k, answer) )
1942 results[numSolnParams++] = answer;
1948 offsets[numNodes] = numSolnParams;
1950 return(FEI_SUCCESS);
1954 int FEDataFilter::getNodalSolution(
int numNodes,
1955 const GlobalID *nodeIDs,
1959 debugOutput(
"FEI: getNodalSolution");
1961 int numActiveNodes = problemStructure_->getNumActiveNodes();
1962 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
1964 if (numActiveNodes <= 0)
return(0);
1966 int numSolnParams = 0;
1972 for(
int i=0; i<numActiveNodes; i++) {
1976 if (offset == numNodes)
break;
1978 GlobalID nodeID = nodeIDs[offset];
1981 offsets[offset++] = numSolnParams;
1988 if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
2002 int numFields = node->getNumFields();
2003 const int* fieldIDs = node->getFieldIDList();
2005 for(
int j=0; j<numFields; j++) {
2006 int size = problemStructure_->getFieldSize(fieldIDs[j]);
2015 for(
int k=0; k<size; k++) {
2016 CHK_ERR( getEqnSolnEntry(thisEqn+k, answer) )
2017 results[numSolnParams++] = answer;
2022 offsets[numNodes] = numSolnParams;
2024 return(FEI_SUCCESS);
2028 int FEDataFilter::getBlockFieldNodeSolution(GlobalID elemBlockID,
2031 const GlobalID *nodeIDs,
2034 debugOutput(
"FEI: getBlockFieldNodeSolution");
2036 int numActiveNodes = problemStructure_->getNumActiveNodes();
2037 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2039 if (numActiveNodes <= 0)
return(0);
2042 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
2044 int fieldSize = problemStructure_->getFieldSize(fieldID);
2045 if (fieldSize <= 0) ERReturn(-1);
2047 if (!block->containsField(fieldID)) {
2048 fei::console_out() <<
"FEDataFilter::getBlockFieldNodeSolution WARNING: fieldID " << fieldID
2049 <<
" not contained in element-block " <<
static_cast<int>(elemBlockID) << FEI_ENDL;
2056 for(
int i=0; i<numNodes; i++) {
2060 GlobalID nodeID = nodeIDs[i];
2067 if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
2083 if (!hasField)
continue;
2085 int offset = fieldSize*i;
2086 for(
int j=0; j<fieldSize; j++) {
2087 double answer = 0.0;
2088 CHK_ERR( getEqnSolnEntry(eqnNumber+j, answer) );
2089 results[offset+j] = answer;
2093 return(FEI_SUCCESS);
2097 int FEDataFilter::getNodalFieldSolution(
int fieldID,
2099 const GlobalID *nodeIDs,
2102 debugOutput(
"FEI: getNodalFieldSolution");
2104 int numActiveNodes = problemStructure_->getNumActiveNodes();
2105 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2107 if (numActiveNodes <= 0)
return(0);
2109 if (problemStructure_->numSlaveEquations() != 0) {
2110 fei::console_out() <<
"FEDataFilter::getEqnSolnEntry ERROR FETI-support is not currently" 2111 <<
" compatible with the FEI's constraint reduction." << FEI_ENDL;
2115 int fieldSize = problemStructure_->getFieldSize(fieldID);
2116 if (fieldSize <= 0) {
2123 for(
int i=0; i<numNodes; i++) {
2127 GlobalID nodeID = nodeIDs[i];
2134 if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
2148 int nodeNumber = node->getNodeNumber();
2155 if (!hasField)
continue;
2157 int dof_id = problemStructure_->getFieldDofMap().get_dof_id(fieldID, 0);
2159 int offset = fieldSize*i;
2161 for(
int j=0; j<fieldSize; j++) {
2162 if (localStartRow_ > eqnNumber || eqnNumber > localEndRow_) {
2163 CHK_ERR( getSharedRemoteSolnEntry(eqnNumber+j, results[offset+j]) );
2167 err = feData_->getSolnEntry(nodeNumber, dof_id+j, results[offset+j]);
2169 fei::console_out() <<
"FEDataFilter::getReducedSolnEntry: nodeNumber " << nodeNumber
2170 <<
" (nodeID " << nodeID <<
"), dof_id "<<dof_id
2171 <<
" couldn't be obtained from FETI on proc " << localRank_ << FEI_ENDL;
2177 return(FEI_SUCCESS);
2181 int FEDataFilter::putBlockNodeSolution(GlobalID elemBlockID,
2183 const GlobalID *nodeIDs,
2185 const double *estimates) {
2187 debugOutput(
"FEI: putBlockNodeSolution");
2189 int numActiveNodes = problemStructure_->getNumActiveNodes();
2191 if (numActiveNodes <= 0)
return(0);
2194 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
2196 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2201 int blk_idx = problemStructure_->getIndexOfBlock(elemBlockID);
2203 for(
int i=0; i<numNodes; i++) {
2205 int err = nodeDB.getNodeWithID(nodeIDs[i], node);
2207 if (err != 0)
continue;
2209 if (!node->hasBlockIndex(blk_idx))
continue;
2211 if (node->getOwnerProc() != localRank_)
continue;
2213 int numFields = node->getNumFields();
2214 const int* fieldIDs = node->getFieldIDList();
2215 const int* fieldEqnNumbers = node->getFieldEqnNumbers();
2217 if (fieldEqnNumbers[0] < localStartRow_ ||
2218 fieldEqnNumbers[0] > localEndRow_)
continue;
2220 int offs = offsets[i];
2222 for(
int j=0; j<numFields; j++) {
2223 int size = problemStructure_->getFieldSize(fieldIDs[j]);
2225 if (block->containsField(fieldIDs[j])) {
2226 for(
int k=0; k<size; k++) {
2229 translateToReducedEqn(fieldEqnNumbers[j]+k, reducedEqn);
2236 return(FEI_SUCCESS);
2240 int FEDataFilter::putBlockFieldNodeSolution(GlobalID elemBlockID,
2243 const GlobalID *nodeIDs,
2244 const double *estimates)
2246 debugOutput(
"FEI: putBlockFieldNodeSolution");
2249 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
2250 if (!block->containsField(fieldID))
return(1);
2252 int fieldSize = problemStructure_->getFieldSize(fieldID);
2253 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2258 std::vector<int> numbers(numNodes);
2263 std::vector<double> data;
2266 if (fieldSize < 1) {
2267 fei::console_out() <<
"FEI Warning, putBlockFieldNodeSolution called for field " 2268 << fieldID<<
", which has size "<<fieldSize<<FEI_ENDL;
2272 numbers.resize(numNodes*fieldSize);
2273 data.resize(numNodes*fieldSize);
2275 catch(std::runtime_error& exc) {
2283 for(
int i=0; i<numNodes; i++) {
2287 if (fieldID < 0) numbers[count++] = node->getNodeNumber();
2291 if (eqn >= localStartRow_ && eqn <= localEndRow_) {
2292 for(
int j=0; j<fieldSize; j++) {
2293 data[count] = estimates[i*fieldSize + j];
2294 problemStructure_->translateToReducedEqn(eqn+j, numbers[count++]);
2302 CHK_ERR( feData_->putNodalFieldData(fieldID, fieldSize,
2303 numNodes, &numbers[0],
2307 return(FEI_SUCCESS);
2311 int FEDataFilter::getBlockElemSolution(GlobalID elemBlockID,
2313 const GlobalID *elemIDs,
2314 int& numElemDOFPerElement,
2321 debugOutput(
"FEI: getBlockElemSolution");
2324 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
2326 std::map<GlobalID,int>& elemIDList = problemStructure_->
2327 getBlockConnectivity(elemBlockID).elemIDs;
2329 int len = block->getNumElements();
2332 if (len > numElems) len = numElems;
2334 numElemDOFPerElement = block->getNumElemDOFPerElement();
2335 std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
2339 if (numElemDOFPerElement <= 0)
return(0);
2341 std::map<GlobalID,int>::const_iterator
2342 elemid_itr = elemIDList.begin();
2344 for(
int i=0; i<len; i++) {
2349 if (elemid_itr->first != elemIDs[i]) {
2350 index = elemid_itr->second;
2353 if (index < 0)
continue;
2355 int offset = i*numElemDOFPerElement;
2357 for(
int j=0; j<numElemDOFPerElement; j++) {
2358 int eqn = elemDOFEqnNumbers[index] + j;
2360 CHK_ERR( getEqnSolnEntry(eqn, answer) )
2362 results[offset+j] = answer;
2368 return(FEI_SUCCESS);
2372 int FEDataFilter::putBlockElemSolution(GlobalID elemBlockID,
2374 const GlobalID *elemIDs,
2376 const double *estimates)
2378 debugOutput(
"FEI: putBlockElemSolution");
2381 CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
2383 std::map<GlobalID,int>& elemIDList = problemStructure_->
2384 getBlockConnectivity(elemBlockID).elemIDs;
2386 int len = block->getNumElements();
2387 if (len > numElems) len = numElems;
2389 int DOFPerElement = block->getNumElemDOFPerElement();
2390 if (DOFPerElement != dofPerElem) {
2391 fei::console_out() <<
"FEI ERROR, putBlockElemSolution called with bad 'dofPerElem' (" 2392 <<dofPerElem<<
"), block "<<elemBlockID<<
" should have dofPerElem==" 2393 <<DOFPerElement<<FEI_ENDL;
2397 std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
2399 if (DOFPerElement <= 0)
return(0);
2401 std::map<GlobalID,int>::const_iterator
2402 elemid_itr = elemIDList.begin();
2404 for(
int i=0; i<len; i++) {
2406 if (elemid_itr->first != elemIDs[i]) {
2407 index = elemid_itr->second;
2410 if (index < 0)
continue;
2412 for(
int j=0; j<DOFPerElement; j++) {
2415 translateToReducedEqn(elemDOFEqnNumbers[i] + j, reducedEqn);
2426 return(FEI_SUCCESS);
2430 int FEDataFilter::getCRMultipliers(
int numCRs,
2432 double* multipliers)
2434 for(
int i=0; i<numCRs; i++) {
2437 multipliers[i] = -999.99;
2444 int FEDataFilter::putCRMultipliers(
int numMultCRs,
2446 const double *multEstimates)
2448 debugOutput(
"FEI: putCRMultipliers");
2450 return(FEI_SUCCESS);
2454 int FEDataFilter::putNodalFieldData(
int fieldID,
2456 const GlobalID* nodeIDs,
2457 const double* nodeData)
2459 debugOutput(
"FEI: putNodalFieldData");
2461 int fieldSize = problemStructure_->getFieldSize(fieldID);
2462 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2464 std::vector<int> nodeNumbers(numNodes);
2466 for(
int i=0; i<numNodes; i++) {
2470 int nodeNumber = node->getNodeNumber();
2471 if (nodeNumber < 0) {
2472 GlobalID nodeID = nodeIDs[i];
2473 fei::console_out() <<
"FEDataFilter::putNodalFieldData ERROR, node with ID " 2474 <<
static_cast<int>(nodeID) <<
" doesn't have an associated nodeNumber " 2475 <<
"assigned. putNodalFieldData shouldn't be called until after the " 2476 <<
"initComplete method has been called." << FEI_ENDL;
2480 nodeNumbers[i] = nodeNumber;
2483 CHK_ERR( feData_->putNodalFieldData(fieldID, fieldSize,
2484 numNodes, &nodeNumbers[0],
2491 int FEDataFilter::putNodalFieldSolution(
int fieldID,
2493 const GlobalID* nodeIDs,
2494 const double* nodeData)
2496 debugOutput(
"FEI: putNodalFieldSolution");
2499 return(putNodalFieldData(fieldID, numNodes, nodeIDs, nodeData));
2502 int fieldSize = problemStructure_->getFieldSize(fieldID);
2503 NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
2505 std::vector<int> eqnNumbers(fieldSize);
2507 for(
int i=0; i<numNodes; i++) {
2513 if (!hasField)
continue;
2521 int FEDataFilter::assembleEqns(
int numRows,
2523 const int* rowNumbers,
2524 const int* colIndices,
2525 const double*
const* coefs,
2526 bool structurallySymmetric,
2529 if (numRows == 0)
return(FEI_SUCCESS);
2531 const int* indPtr = colIndices;
2532 for(
int i=0; i<numRows; i++) {
2533 int row = rowNumbers[i];
2535 const double* coefPtr = coefs[i];
2537 CHK_ERR(giveToMatrix(1, &row, numCols, indPtr, &coefPtr, mode));
2539 if (!structurallySymmetric) indPtr += numCols;
2542 return(FEI_SUCCESS);
2546 int FEDataFilter::assembleRHS(
int numValues,
2548 const double* coefs,
2554 if (problemStructure_->numSlaveEquations() == 0) {
2555 CHK_ERR( giveToRHS(numValues, coefs, indices, mode) );
2556 return(FEI_SUCCESS);
2559 for(
int i = 0; i < numValues; i++) {
2560 int eqn = indices[i];
2561 if (eqn < 0)
continue;
2563 CHK_ERR( giveToRHS(1, &(coefs[i]), &eqn, mode ) );
2566 return(FEI_SUCCESS);
2570 void FEDataFilter::debugOutput(
const char* mesg)
2572 if (Filter::logStream() != NULL) {
2573 (*logStream()) << mesg << FEI_ENDL;
int getParameters(int &numParams, char **¶mStrings)
bool getFieldEqnNumber(int fieldID, int &eqnNumber) const
std::vector< int > & getMasterFieldIDs()
virtual void parameters(const fei::ParameterSet ¶mset)
void getFieldID(int eqnNumber, int &fieldID, int &offset_into_field) const
int getNumNodeDescriptors() const
int getNodeWithID(GlobalID nodeID, const NodeDescriptor *&node) const
void getNodeAtIndex(int i, const NodeDescriptor *&node) const
int ** fieldIDsTablePtr()
int binarySearch(const T &item, const T *list, int len)
std::ostream & console_out()
int localProc(MPI_Comm comm)
int getNodeWithNumber(int nodeNumber, const NodeDescriptor *&node) const
int getNodeWithEqn(int eqnNumber, const NodeDescriptor *&node) const
std::vector< int > & getMasters()
int numProcs(MPI_Comm comm)