diff --git a/Solver/dofManager.h b/Solver/dofManager.h index 9a1ae18f770a929693e28321e8597cef4b05d508..778259f889a8a75908624845aa973db8d28b2037 100644 --- a/Solver/dofManager.h +++ b/Solver/dofManager.h @@ -11,6 +11,8 @@ #include <map> #include "MVertex.h" #include "linearSystem.h" +#include "fullMatrix.h" + class Dof{ private: @@ -18,6 +20,7 @@ class Dof{ long int _entity; // "i": node, edge, group, etc. int _type; // "f": basis function type index, etc. public: + Dof(void) {} // eric Dof(long int entity, int type) : _entity(entity), _type(type) {} inline long int getEntity() const { return _entity; } inline int getType() const { return _type; } @@ -133,6 +136,43 @@ class dofManager{ } } } + + inline void assemble(const std::vector<Dof> &dofsR ,const std::vector<Dof> &dofsC, const fullMatrix<dataMat> &localMatrix) + { + if (!_current->isAllocated()) _current->allocate(unknown.size()); + const int nbR = localMatrix.size1(); + const int nbC = localMatrix.size2(); + std::vector<int> tabR(nbR); + std::vector<int> tabC(nbC); + for (int R=0;R<nbR;++R) + { + std::map<Dof, int>::iterator it = unknown.find(dofsR[R]); + if (it != unknown.end()) tabR[R]=it->second; else tabR[R]=-1; + } + for (int C=0;C<nbC;++C) + { + std::map<Dof, int>::iterator it = unknown.find(dofsC[C]); + if (it != unknown.end()) tabC[C]=it->second; else tabC[C]=-1; + } + + for (int R=0;R<nbR;++R) { + for (int C=0;C<nbC;++C) { + if (tabR[R] != -1) { + if (tabC[C] !=-1) { + _current->addToMatrix(tabR[R], tabC[C], localMatrix(R,C)); + } + else { + typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(dofsC[C]); + if (itFixed != fixed.end()) { + _current->addToRightHandSide(R, -localMatrix(R,C) * itFixed->second); + } + } + } + } + } + } + + inline void assemble(int entR, int typeR, int entC, int typeC, const dataMat &value) { assemble(Dof(entR, typeR), Dof(entC, typeC), value); diff --git a/Solver/femTerm.h b/Solver/femTerm.h index 65c6938f5fa719d010b1c409b994ea7403ab7d92..ff12b9dcb1079e0a3b4bc12b931dce6841ece6c9 100644 --- a/Solver/femTerm.h +++ b/Solver/femTerm.h @@ -75,13 +75,24 @@ class femTerm { { const int nbR = localMatrix.size1(); const int nbC = localMatrix.size2(); - for (int j = 0; j < nbR; j++){ +/* + for (int j = 0; j < nbR; j++) + { Dof R = getLocalDofR(se, j); - for (int k = 0; k < nbC; k++){ + for (int k = 0; k < nbC; k++) + { Dof C = getLocalDofC(se, k); dm.assemble(R, C, localMatrix(j, k)); } } +*/ + std::vector<Dof> tabR(nbR); + std::vector<Dof> tabC(nbC); + for (int R = 0; R < nbR; ++R) + tabR[R]=getLocalDofR(se, R); + for (int C = 0; C < nbC; ++C) + tabC[C]=getLocalDofC(se, C); + dm.assemble(tabR,tabC,localMatrix); } void dirichletNodalBC(int physical, int dim, int comp, int field, const simpleFunction<dataVec> &e,