diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt index e88ccd212ef704586aa5b5495ecce99a85fb8921..83714dfa9e466ad9ff5c1c7f23aebe4830ab007f 100644 --- a/Solver/CMakeLists.txt +++ b/Solver/CMakeLists.txt @@ -13,6 +13,7 @@ set(SRC multiscaleLaplace.cpp dgGroupOfElements.cpp dgAlgorithm.cpp + dgMgAlgorithm.cpp dgRungeKutta.cpp dgRungeKuttaMultirate.cpp dgConservationLaw.cpp diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp index cbf6afa0af39791c8142f6544b84ae34ac80c3d6..32df24e05fe232ab9d5e454ef4508f4cf8b2cc0b 100644 --- a/Solver/dgDofContainer.cpp +++ b/Solver/dgDofContainer.cpp @@ -242,6 +242,44 @@ void dgDofContainer::L2Projection(std::string functionName){ } } +void dgDofContainer::Mesh2Mesh_BuildL2Projection(linearSystemCSR<double> &projector,dgDofContainer &donor){ +// projector.allocate(); +// projector.addToMatrix(,,); +/*nt nbCoarseNodes = 0; // OLD CODE + int nbCoarseElements = 0; + int r = 0; + for (int kFine = 0;kFine < fineGroups.getNbElementGroups();kFine++) { + int nbFineNodes = fineGroups.getElementGroup(kFine)->getNbNodes(); + fineSol.getGroupProxy(kFine).scale(0); + for (int fineElement = 0 ; fineElement<fineGroups.getElementGroup(kFine)->getNbElements() ;++fineElement) { + for (int fineNodes = 0 ; fineNodes<nbFineNodes ;++fineNodes) { + int kCoarse = 0; + int c0 = 0; + for (int i = startIndex[r++]; i < startIndex[r] - 1; i++){ + int c = columns[i]; + for (;c-c0 < nbCoarseNodes*nbCoarseElements;kCoarse++) { + c0 += nbCoarseNodes*nbCoarseElements; + nbCoarseNodes = coarseGroups.getElementGroup(kCoarse)->getNbNodes(); + nbCoarseElements = coarseGroups.getElementGroup(kCoarse)->getNbElements(); + } + int coarseElement = (c-c0)/nbCoarseNodes; + int coarseNodes = c-c0-nbCoarseNodes*coarseElement; + for (int m = 0; m < nbFields; m++){ + (fineSol.getGroupProxy(kFine))(fineNodes,fineElement*nbFields+m) += values[i]*(coarseSol.getGroupProxy(kCoarse))(coarseNodes,coarseElement*nbFields+m); + } + } + } + }*/ +} + +void dgDofContainer::Mesh2Mesh_ApplyL2Projection(linearSystemCSR<double> &projector,dgDofContainer &donor){ + scale(0.); + int *startIndex; + int *columns; + double *values; + projector.getMatrix(startIndex,columns,values); +} + void dgDofContainer::exportGroupIdMsh(){ // the elementnodedata::export does not work !! diff --git a/Solver/dgDofContainer.h b/Solver/dgDofContainer.h index ad926cd1ef44d70dad7215c00b556a5552b066f7..29318b62f1973f7761159ff062f3b0ba313edaae 100644 --- a/Solver/dgDofContainer.h +++ b/Solver/dgDofContainer.h @@ -3,6 +3,7 @@ #include <vector> #include "fullMatrix.h" +#include "linearSystemCSR.h" class dgGroupCollection; class dgGroupOfElements; @@ -42,7 +43,8 @@ public: void load(const std::string name); void setAll(double v); void L2Projection(std::string functionName); - void Mesh2Mesh_L2Projection(dgDofContainer &other); + void Mesh2Mesh_BuildL2Projection(linearSystemCSR<double> &projector,dgDofContainer &donor); + void Mesh2Mesh_ApplyL2Projection(linearSystemCSR<double> &projector,dgDofContainer &donor); void exportMsh(const std::string name); void exportGroupIdMsh(); void exportMultirateGroupMsh(); diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp index 7248b24dd5d09875e59f15f8869fa708f858a2a9..192bc6a794915608df899d848f41c4668d5718e0 100644 --- a/Solver/dgGroupOfElements.cpp +++ b/Solver/dgGroupOfElements.cpp @@ -168,6 +168,7 @@ dgGroupOfElements::~dgGroupOfElements(){ delete _elementVolume; } + dgGroupOfFaces::~dgGroupOfFaces() { if (!_faces.size())return; diff --git a/Solver/function.cpp b/Solver/function.cpp index 034d02560952d866d697aa0fa1f4e24b3532400f..1bd3041b8a311f3982eea8f1066a64e8fdd0e58e 100644 --- a/Solver/function.cpp +++ b/Solver/function.cpp @@ -291,7 +291,7 @@ public: for (int k=0;k<_dofContainer->getNbFields();k++){ _value(i,k) = 0.0; for (int j=0;j<fSize;j++){ - _value(i,k) += solEl(j,k)*fs[j]; + _value(i,k) += solEl(j,k)*fs[j]; } } } diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp index fd31746723c5fbe2ddd3beffb6af44b6855a8925..570657b5c450bc6e952f7027f9b458f6f8c81b86 100644 --- a/Solver/linearSystemCSR.cpp +++ b/Solver/linearSystemCSR.cpp @@ -282,6 +282,31 @@ void sortColumns_(int NbLines, delete[] count; } +template<> + void linearSystemCSR<double>::getMatrix(INDEX_TYPE*& jptr,INDEX_TYPE*& ai,double*& a) + { + jptr = (INDEX_TYPE*) _jptr->array; + ai = (INDEX_TYPE*) _ai->array; + a = ( double * ) _a->array; + if (!sorted) + sortColumns_(_b->size(), + CSRList_Nbr(_a), + (INDEX_TYPE *) _ptr->array, + jptr, + ai, + a); + sorted = true; + } + +#include "Bindings.h" + +template<> + void linearSystemCSR<double>::registerBindings(binding *b) + { + classBinding *cb = b->addClass< linearSystemCSR<double> >("linearSystemCSRdouble"); +// cb->setDescription("A GModel contains a geometrycal and it's mesh."); + } + #if defined(HAVE_GMM) #include "gmm.h" diff --git a/Solver/linearSystemCSR.h b/Solver/linearSystemCSR.h index 84e4c13e4c2e9da4072255c2c742f3c8b9776fd4..58e343741af4e59b276d249d98f5fb3c815caa37 100644 --- a/Solver/linearSystemCSR.h +++ b/Solver/linearSystemCSR.h @@ -11,6 +11,8 @@ #include "GmshMessage.h" #include "linearSystem.h" +class binding; + typedef int INDEX_TYPE ; typedef struct { int nmax; @@ -83,6 +85,8 @@ class linearSystemCSR : public linearSystem<scalar> { } else ptr[position] = n; } + virtual void getMatrix(INDEX_TYPE*& jptr,INDEX_TYPE*& ai,double*& a); + virtual scalar getFromMatrix (int row, int col) const { Msg::Error("getFromMatrix not implemented for CSR"); @@ -110,6 +114,7 @@ class linearSystemCSR : public linearSystem<scalar> { { for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.; } + static void registerBindings(binding *b); }; template <class scalar>