diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 667d49dae25df24f22b2eddf16a80e64686eeaed..809a3885bfd25b1d652d45304765273a4ddb4a49 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -520,7 +520,7 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
MElement *copy = factory.create(eType, vmv, ++numEle, ePart);
double **nodeLs = new double*[e->getNumPrimaryVertices()];
- for(int i = 0; i < e->getNumPrimaryVertices(); i++) nodeLs[i] = new double[verticesLs.size2()];
+ // for(int i = 0; i < e->getNumPrimaryVertices(); i++) nodeLs[i] = new double[verticesLs.size2()];
switch (eType) {
case MSH_TET_4 :
@@ -534,7 +534,10 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z());
- for(int i = 0; i < 4; i++) nodeLs[i] = &verticesLs(e->getVertex(i)->getIndex(),0);
+ for(int i = 0; i < 4; i++) {
+ nodeLs[i] = &verticesLs(e->getVertex(i)->getIndex(),0);
+ //for(int j = 0; j < verticesLs.size2(); j++)nodeLs(i,j) = verticesLs[j];
+ }
isCut = T.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
tetras, quads, triangles, 0, nodeLs);
}
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 18d0efef5faedf1ff334c6e59802924f13d38157..4fe957c33d55f753f0c4718642a56924679b3d8c 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -135,6 +135,12 @@ class fullMatrix
}
int gemm(lua_State *L);
#endif // HAVE LUA
+ fullMatrix(scalar *original, int r, int c){
+ _r = r;
+ _c = c;
+ _own_data = false;
+ _data = original;
+ }
fullMatrix(fullMatrix<scalar> &original, int c_start, int c){
_c = c;
_r = original._r;
diff --git a/Post/PView.cpp b/Post/PView.cpp
index 8f59d995badd2c3301e26709534c7abc7cc5ab4f..7595047c5b6f23dc280320ce18aaddf24c44fb30 100644
--- a/Post/PView.cpp
+++ b/Post/PView.cpp
@@ -112,7 +112,7 @@ PView::PView(std::string xname, std::string yname,
PView::PView(std::string name, std::string type,
GModel *model, std::map<int, std::vector<double> > &data,
- double time)
+ double time, int numC)
{
_init();
PViewDataGModel::DataType t;
@@ -127,7 +127,7 @@ PView::PView(std::string name, std::string type,
return;
}
PViewDataGModel *d = new PViewDataGModel(t);
- d->addData(model, data, 0, time, 1);
+ d->addData(model, data, 0, time, 1, numC);
d->setName(name);
d->setFileName(name + ".msh");
_data = d;
@@ -138,10 +138,10 @@ PView::PView(std::string name, std::string type,
}
void PView::addStep(GModel *model, std::map<int, std::vector<double> > &data,
- double time)
+ double time, int numC)
{
PViewDataGModel *d = dynamic_cast<PViewDataGModel*>(_data);
- if(d) d->addData(model, data, d->getNumTimeSteps(), time, 1);
+ if(d) d->addData(model, data, d->getNumTimeSteps(), time, 1, numC);
else Msg::Error("Can only add step data to model-based datasets");
}
diff --git a/Post/PView.h b/Post/PView.h
index d621a93e4e1b6d6ddc815739da0f303db652a45d..04d8727acd6953072d10d6979d58ddc6274012d9 100644
--- a/Post/PView.h
+++ b/Post/PView.h
@@ -52,10 +52,11 @@ class PView{
std::vector<double> &x, std::vector<double> &y);
// construct a new model-based view from a bunch of data
PView(std::string name, std::string type, GModel *model,
- std::map<int, std::vector<double> > &data, double time=0.);
+ std::map<int, std::vector<double> > &data, double time=0.,
+ int numComp = -1);
// add a new time step to a given model-based view
void addStep(GModel *model, std::map<int, std::vector<double> > &data,
- double time=0.);
+ double time=0.,int numComp = -1);
// default destructor
~PView();
diff --git a/Post/PViewDataGModel.h b/Post/PViewDataGModel.h
index 9a5cf29edf2b58f68e090f9c8de8f7a10a2221b1..2cc0c60fc3410d0b385241bc3e0703049d7d0d37 100644
--- a/Post/PViewDataGModel.h
+++ b/Post/PViewDataGModel.h
@@ -185,7 +185,7 @@ class PViewDataGModel : public PViewData {
// node or element number depending on the type of dataset; all the
// vectors are supposed to have the same length)
bool addData(GModel *model, std::map<int, std::vector<double> > &data,
- int step, double time, int partition);
+ int step, double time, int partition, int numC);
// I/O routines
bool readMSH(std::string fileName, int fileIndex, FILE *fp, bool binary,
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index 6979d617c9e18e450c589afcde5a04f170902f8e..60bfa2c41ec216e1a33e4d2aa19b63dbad8c59c0 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -12,14 +12,17 @@
#include "StringUtils.h"
bool PViewDataGModel::addData(GModel *model, std::map<int, std::vector<double> > &data,
- int step, double time, int partition)
+ int step, double time, int partition, int numC)
{
if(data.empty()) return false;
int numComp = 9;
- for(std::map<int, std::vector<double> >::iterator it = data.begin();
- it != data.end(); it++)
- numComp = std::min(numComp, (int)it->second.size());
+ if (numC < 0){
+ for(std::map<int, std::vector<double> >::iterator it = data.begin();
+ it != data.end(); it++)
+ numComp = std::min(numComp, (int)it->second.size());
+ }
+ else numComp = numC;
while(step >= (int)_steps.size())
_steps.push_back(new stepData<double>(model, numComp));
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index 60898af6e31f95eb5e2e0deb4f1e2c8e654ea68a..b230ce42754bb5722167e645c96d87eb117c0c5f 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -20,6 +20,7 @@ set(SRC
dgConservationLawWaveEquation.cpp
function.cpp
functionDerivator.cpp
+ luaFunction.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Solver *.h)
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 9d24b827039ae1a8c7829590ab51b6e9d88251b3..0203bcf63ad39d616755ab5258b3600845fdc2fd 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -17,7 +17,7 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
const dgGroupOfElements & group,
const fullMatrix<double> &solution,
fullMatrix<double> &residual // the residual
- )
+ )
{
// ----- 1 ---- get the solution at quadrature points
// ----- 1.1 --- allocate a matrix of size (nbFields * nbElements, nbQuadraturePoints)
@@ -114,7 +114,8 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
}
if(sourceTerm){
residual.gemm(group.getSourceRedistributionMatrix(),Source);
- delete sourceTerm;
+ // FIXME (JF) : for now TEST TEST
+ // delete sourceTerm;
}
}
@@ -232,49 +233,52 @@ void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
}
}
- void dgAlgorithm::rungeKutta (
- const dgConservationLaw &claw, // conservation law
- std::vector<dgGroupOfElements*> &eGroups, // group of elements
- std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
- std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
- double h, // time-step
- fullMatrix<double> &residu,
- fullMatrix<double> &sol,
- int orderRK) // order of RK integrator
-{
+
+ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation law
+ std::vector<dgGroupOfElements*> &eGroups, // group of elements
+ std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
+ std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
+ double h, // time-step
+ std::vector<fullMatrix<double> *> &residu,
+ std::vector<fullMatrix<double> *> &sol,
+ int orderRK) // order of RK integrator
+ {
// U_{n+1}=U_n+h*(SUM_i a_i*K_i)
// K_i=M^(-1)R(U_n+b_i*K_{i-1})
- double a[4] = {h/6.0,h/3.0,h/3.0,h/6.0};
- double b[4] = {0.,h/2.0,h/2.0,h};
-
- fullMatrix<double> K(sol);
- // Current updated solution
- fullMatrix<double> Unp(sol);
- fullMatrix<double> residuEl, KEl;
- fullMatrix<double> iMassEl;
-
- int nbNodes = eGroups[0]->getNbNodes();
- int nbFields = sol.size2()/eGroups[0]->getNbElements();
-
- for(int j=0; j<orderRK;j++){
- if(j!=0){
- K.scale(b[j]);
- K.add(sol);
- }
- this->residual(claw,eGroups,fGroups,bGroups,K,residu);
- K.scale(0.);
- for(int i=0;i<eGroups[0]->getNbElements();i++) {
- residuEl.setAsProxy(residu,i*nbFields,nbFields);
- KEl.setAsProxy(K,i*nbFields,nbFields);
- iMassEl.setAsProxy(eGroups[0]->getInverseMassMatrix(),i*nbNodes,nbNodes);
- iMassEl.mult(residuEl,KEl);
- }
- Unp.add(K,a[j]);
- }
- sol=Unp;
-}
+ /*
+ double a[4] = {h/6.0,h/3.0,h/3.0,h/6.0};
+ double b[4] = {0.,h/2.0,h/2.0,h};
+
+ fullMatrix<double> K(sol);
+ // Current updated solution
+ fullMatrix<double> Unp(sol);
+ fullMatrix<double> residuEl, KEl;
+ fullMatrix<double> iMassEl;
+
+ int nbNodes = eGroups[0]->getNbNodes();
+ int nbFields = sol.size2()/eGroups[0]->getNbElements();
+
+ for(int j=0; j<orderRK;j++){
+ if(j!=0){
+ K.scale(b[j]);
+ K.add(sol);
+ }
+ this->residual(claw,eGroups,fGroups,bGroups,K,residu);
+ K.scale(0.);
+ for(int i=0;i<eGroups[0]->getNbElements();i++) {
+ residuEl.setAsProxy(residu,i*nbFields,nbFields);
+ KEl.setAsProxy(K,i*nbFields,nbFields);
+ iMassEl.setAsProxy(eGroups[0]->getInverseMassMatrix(),i*nbNodes,nbNodes);
+ iMassEl.mult(residuEl,KEl);
+ }
+ Unp.add(K,a[j]);
+ }
+ sol=Unp;
+ }
+ */
+ }
void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
@@ -386,37 +390,47 @@ void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
}
}
-// works only if there is only 1 group of element
+// works for any number of groups
void dgAlgorithm::residual( const dgConservationLaw &claw,
- std::vector<dgGroupOfElements*> &eGroups, //group of elements
- std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
- std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
- fullMatrix<double> &solution, // solution
- fullMatrix<double> &residu) // residual
+ std::vector<dgGroupOfElements*> &eGroups, //group of elements
+ std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
+ std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
+ std::vector<fullMatrix<double> *> &solution, // solution
+ std::vector<fullMatrix<double> *> &residu) // residual
{
int nbFields=claw.nbFields();
- dgAlgorithm algo;
- residu.scale(0);
//volume term
for(size_t i=0;i<eGroups.size() ; i++) {
- algo.residualVolume(claw,*eGroups[i],solution,residu);
+ residu[i]->scale(0);
+ residualVolume(claw,*eGroups[i],*solution[i],*residu[i]);
}
//interface term
for(size_t i=0;i<fGroups.size() ; i++) {
dgGroupOfFaces &faces = *fGroups[i];
+ int iGroupLeft = -1, iGroupRight = -1;
+ for(size_t j=0;j<eGroups.size() ; j++) {
+ if (eGroups[i] == &faces.getGroupLeft())iGroupLeft = j;
+ if (eGroups[i] == &faces.getGroupRight())iGroupRight= j;
+ }
+
fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
- faces.mapToInterface(nbFields, solution, solution, solInterface);
- algo.residualInterface(claw,faces,solInterface,solution,solution,residuInterface);
- faces.mapFromInterface(nbFields, residuInterface, residu, residu);
+ faces.mapToInterface(nbFields, *solution[iGroupLeft], *solution[iGroupRight], solInterface);
+ residualInterface(claw,faces,solInterface,*solution[iGroupLeft], *solution[iGroupRight],residuInterface);
+ faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupLeft]);
}
//boundaries
for(size_t i=0;i<bGroups.size() ; i++) {
dgGroupOfFaces &faces = *bGroups[i];
+ int iGroupLeft = -1, iGroupRight = -1;
+ for(size_t j=0;j<eGroups.size() ; j++) {
+ if (eGroups[i] == &faces.getGroupLeft())iGroupLeft = j;
+ if (eGroups[i] == &faces.getGroupRight())iGroupRight= j;
+ }
fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
- faces.mapToInterface(nbFields, solution, solution, solInterface);
- algo.residualBoundary(claw,faces,solInterface,solution,residuInterface);
- faces.mapFromInterface(nbFields, residuInterface, residu, residu);
+ faces.mapToInterface(nbFields, *solution[iGroupLeft], *solution[iGroupRight], solInterface);
+ residualBoundary(claw,faces,solInterface,*solution[iGroupLeft],residuInterface);
+ faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupRight]);
}
}
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index f51f050e423634c154b1258b472f86b835b0a0ef..dfad3132e488ba949dd3a9674a8e4c4e08e69ca1 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -11,50 +11,50 @@ class dgTerm;
class dgAlgorithm {
public :
- static void residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
- const dgConservationLaw &claw, // the conservation law
- const dgGroupOfElements & group,
- const fullMatrix<double> &solution,
- fullMatrix<double> &residual);
- static void residualInterface ( //dofManager &dof, // the DOF manager (maybe useless here)
- const dgConservationLaw &claw, // the conservation law
- const dgGroupOfFaces & group,
- const fullMatrix<double> &solution, // ! at face nodes
- fullMatrix<double> &solutionRight,
- fullMatrix<double> &solutionLeft,
- fullMatrix<double> &residual); // ! at face nodes
- static void residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
- const dgConservationLaw &claw, // the conservation law
- const dgGroupOfFaces &group,
- const fullMatrix<double> &solution, // solution !! at faces nodes
- fullMatrix<double> &solutionLeft,
- fullMatrix<double> &residual // residual !! at faces nodes
- );
- // works only if there is only 1 group of element
+ static void residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
+ const dgConservationLaw &claw, // the conservation law
+ const dgGroupOfElements & group,
+ const fullMatrix<double> &solution,
+ fullMatrix<double> &residual);
+ static void residualInterface ( //dofManager &dof, // the DOF manager (maybe useless here)
+ const dgConservationLaw &claw, // the conservation law
+ const dgGroupOfFaces & group,
+ const fullMatrix<double> &solution, // ! at face nodes
+ fullMatrix<double> &solutionRight,
+ fullMatrix<double> &solutionLeft,
+ fullMatrix<double> &residual); // ! at face nodes
+ static void residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
+ const dgConservationLaw &claw, // the conservation law
+ const dgGroupOfFaces &group,
+ const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &solutionLeft,
+ fullMatrix<double> &residual // residual !! at faces nodes
+ );
+ // works only if there is only 1 group of element
static void residual( const dgConservationLaw &claw,
- std::vector<dgGroupOfElements*> &eGroups, //group of elements
- std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
- std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
- fullMatrix<double> &solution, // solution !! at faces nodes
- fullMatrix<double> &residual // residual !! at faces nodes
- );
+ std::vector<dgGroupOfElements*> &eGroups, //group of elements
+ std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
+ std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
+ std::vector<fullMatrix<double> *> &solution, // solution !! at faces nodes
+ std::vector<fullMatrix<double> *> &residual // residual !! at faces nodes
+ );
void rungeKutta (
- const dgConservationLaw &claw,
- std::vector<dgGroupOfElements*> &eGroups, //group of elements
- std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
- std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
- double h,
- fullMatrix<double> &residu,
- fullMatrix<double> &sol,
- int orderRK=4);
+ const dgConservationLaw &claw,
+ std::vector<dgGroupOfElements*> &eGroups, //group of elements
+ std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
+ std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
+ double h,
+ std::vector<fullMatrix<double> *> &residu,
+ std::vector<fullMatrix<double> *> &sol,
+ int orderRK=4);
static void multAddInverseMassMatrix ( /*dofManager &dof,*/
- const dgGroupOfElements & group,
- fullMatrix<double> &residu,
- fullMatrix<double> &sol);
+ const dgGroupOfElements & group,
+ fullMatrix<double> &residu,
+ fullMatrix<double> &sol);
static void buildGroups(GModel *model, int dim, int order,
- std::vector<dgGroupOfElements*> &eGroups,
- std::vector<dgGroupOfFaces*> &fGroups,
- std::vector<dgGroupOfFaces*> &bGroups);
+ std::vector<dgGroupOfElements*> &eGroups,
+ std::vector<dgGroupOfFaces*> &fGroups,
+ std::vector<dgGroupOfFaces*> &bGroups);
};
diff --git a/Solver/dgMainLua.cpp b/Solver/dgMainLua.cpp
index 487487f38c36ab07a62f64c1d42aa6f84815c33e..cb0ac1845b3172ecdb43fa79ba7a119aee64dbe2 100644
--- a/Solver/dgMainLua.cpp
+++ b/Solver/dgMainLua.cpp
@@ -2,8 +2,8 @@
#include "Gmsh.h"
#include "LuaBindings_Geo.h"
#include "dgSystemOfEquations.h"
+#include "luaFunction.h"
#include "function.h"
-
void report_errors(lua_State *L, int status)
{
if ( status!=0 ) {
@@ -28,6 +28,7 @@ int main(int argc, char *argv[])
dgSystemOfEquations::Register(L);
fullMatrix<double>::Register(L);
function::registerDefaultFunctions();
+ RegisterFunctions(L);
int s = luaL_loadfile(L, argv[1]);
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 4faccae24e068cf76a2d749249e5551596abd1fe..6cc5af311124f95cf17de611f4f7a66a5cb24db8 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -1,5 +1,11 @@
+#include <stdio.h>
+#include <stdlib.h>
#include "dgSystemOfEquations.h"
#include "LuaBindings_Geo.h"
+#include "function.h"
+#include "MElement.h"
+#include "PView.h"
+#include "PViewData.h"
#if defined(HAVE_LUA)
// define the name of the object that will be used in Lua
@@ -15,6 +21,7 @@ Luna<dgSystemOfEquations>::RegType dgSystemOfEquations::methods[] = {
_method(dgSystemOfEquations, setOrder),
_method(dgSystemOfEquations, setup),
_method(dgSystemOfEquations, addBoundaryCondition),
+ _method(dgSystemOfEquations, L2Projection),
{0,0}
};
@@ -24,6 +31,22 @@ void dgSystemOfEquations::Register (lua_State *L){
Luna<dgSystemOfEquations>::Register(L);
}
+class dgConservationLawL2Projection : public dgConservationLaw {
+ std::string _functionName;
+public:
+ dgConservationLawL2Projection(const std::string & functionName,
+ dgConservationLaw &_claw) :
+ _functionName(functionName)
+ {
+ _nbf =_claw.nbFields();
+ }
+ dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const {
+ //return new gaussian(cacheMap,0.2,0.3);
+ return &cacheMap.get(_functionName);
+ }
+};
+
+
// system of equations is build using a mesh
dgSystemOfEquations::dgSystemOfEquations(lua_State *L){
// get the number of arguments
@@ -71,6 +94,14 @@ int dgSystemOfEquations::addBoundaryCondition (lua_State *L){
else throw;
}
+// do a L2 projection
+int dgSystemOfEquations::L2Projection (lua_State *L){
+ dgConservationLawL2Projection Law(std::string(luaL_checkstring(L, 1)),*_claw);
+ for (int i=0;i<_elementGroups.size();i++){
+ _algo->residualVolume(Law,*_elementGroups[i],*_solutionProxys[i],*_rightHandSideProxys[i]);
+ _algo->multAddInverseMassMatrix(*_elementGroups[i],*_rightHandSideProxys[i],*_solutionProxys[i]);
+ }
+}
// ok, we can setup the groups and create solution vectors
int dgSystemOfEquations::setup(lua_State *L){
@@ -81,13 +112,28 @@ int dgSystemOfEquations::setup(lua_State *L){
_elementGroups,
_faceGroups,
_boundaryGroups);
- //for now, we suppose there is only one group of elements
- int nbNodes = _elementGroups[0]->getNbNodes();
- int nbElements = _elementGroups[0]->getNbElements();
- // allocate solution and RHS
+ // compute the total size of the soution
+ _dataSize = 0;
int nbFields = _claw->nbFields();
- _solution = new fullMatrix<double> (nbNodes,nbFields*nbElements);
- _rightHandSide = new fullMatrix<double> (nbNodes,nbFields*nbElements);
+ for (int i=0;i<_elementGroups.size();i++){
+ int nbNodes = _elementGroups[i]->getNbNodes();
+ int nbElements = _elementGroups[i]->getNbElements();
+ _dataSize += nbNodes*nbFields*nbElements;
+ }
+ // allocate the big vectors
+ _solution = new double [_dataSize];
+ _rightHandSide = new double [_dataSize];
+ // create proxys for each group
+ int offset = 0;
+ _solutionProxys.resize(_elementGroups.size());
+ _rightHandSideProxys.resize(_elementGroups.size());
+ for (int i=0;i<_elementGroups.size();i++){
+ int nbNodes = _elementGroups[i]->getNbNodes();
+ int nbElements = _elementGroups[i]->getNbElements();
+ _solutionProxys[i] = new fullMatrix<double> (&_solution[offset*sizeof(double)],nbNodes, nbFields*nbElements);
+ _rightHandSideProxys[i] = new fullMatrix<double> (&_rightHandSide[offset*sizeof(double)],nbNodes, nbFields*nbElements);
+ offset += nbNodes*nbFields*nbElements;
+ }
}
@@ -97,7 +143,7 @@ int dgSystemOfEquations::getSolution(lua_State *L){
}
int dgSystemOfEquations::computeRHS(lua_State *L){
- _algo->residual(*_claw, _elementGroups, _faceGroups, _boundaryGroups, *_solution, *_rightHandSide);
+ _algo->residual(*_claw, _elementGroups, _faceGroups, _boundaryGroups, _solutionProxys, _rightHandSideProxys);
lua_pushlightuserdata (L, _solution);
return 1;
}
@@ -108,35 +154,78 @@ int dgSystemOfEquations::getRHS(lua_State *L){
}
int dgSystemOfEquations::exportSolution(lua_State *L){
+ std::string outputFile(luaL_checkstring(L, 1));
+ export_solution_as_is(outputFile);
}
#endif // HAVE_LUA
-dgSystemOfEquations::dgSystemOfEquations(GModel *obj, dgConservationLaw *claw,int order){
- setup(obj,claw,order);
-}
-
dgSystemOfEquations::~dgSystemOfEquations(){
- delete _solution;
- delete _rightHandSide;
+ for (int i=0;i<_elementGroups.size();i++){
+ delete _solutionProxys[i];
+ delete _rightHandSideProxys[i];
+ delete _elementGroups[i];
+ }
+ delete [] _solution;
+ delete [] _rightHandSide;
}
-void dgSystemOfEquations::setup(GModel *obj, dgConservationLaw *claw, int order){
- _gm = obj;
- _claw = claw;
- _algo = new dgAlgorithm;
- _order= order;
- _dimension = _gm->getNumRegions() ? 3 : 2;
- _algo->buildGroups(_gm,
- _dimension,
- _order,
- _elementGroups,
- _faceGroups,
- _boundaryGroups);
- //for now, we suppose there is only one group of elements
- int nbNodes = _elementGroups[0]->getNbNodes();
- int nbElements = _elementGroups[0]->getNbElements();
- // allocate solution and RHS
- int nbFields = _claw->nbFields();
- _solution = new fullMatrix<double> (nbNodes,nbFields*nbElements);
- _rightHandSide = new fullMatrix<double> (nbNodes,nbFields*nbElements);
+void dgSystemOfEquations::export_solution_as_is (const std::string &name){
+ // the elementnodedata::export does not work !!
+
+ for (int ICOMP = 0; ICOMP<_claw->nbFields();++ICOMP){
+ char aaa[245];
+ sprintf(aaa,"%s_COMP_%d.msh",name.c_str(),ICOMP);
+ FILE *f = fopen (aaa,"w");
+ int COUNT = 0;
+ for (int i=0;i < _elementGroups.size() ;i++){
+ COUNT += _elementGroups[i]->getNbElements();
+ }
+ fprintf(f,"$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
+ fprintf(f,"$ElementNodeData\n");
+ fprintf(f,"1\n");
+ fprintf(f,"\"%s\"\n",name.c_str());
+ fprintf(f,"1\n");
+ fprintf(f,"0.0\n");
+ fprintf(f,"3\n");
+ fprintf(f,"0\n 1\n %d\n",COUNT);
+ for (int i=0;i < _elementGroups.size() ;i++){
+ for (int iElement=0 ; iElement<_elementGroups[i]->getNbElements() ;++iElement) {
+ MElement *e = _elementGroups[i]->getElement(iElement);
+ int num = e->getNum();
+ fullMatrix<double> sol = getSolutionProxy (i, iElement);
+ fprintf(f,"%d %d",num,sol.size1());
+ for (int k=0;k<sol.size1();++k) {
+ fprintf(f,"%12.5E ",sol(k,ICOMP));
+ }
+ fprintf(f,"\n");
+ }
+ }
+ fprintf(f,"$EndElementNodeData\n");
+ fclose(f);
+ }
+ return;
+ // we should discuss that : we do a copy of the solution, this should
+ // be avoided !
+
+ std::map<int, std::vector<double> > data;
+
+ for (int i=0;i < _elementGroups.size() ;i++){
+ for (int iElement=0 ; iElement<_elementGroups[i]->getNbElements() ;++iElement) {
+ MElement *e = _elementGroups[i]->getElement(iElement);
+ int num = e->getNum();
+ fullMatrix<double> sol = getSolutionProxy (i, iElement);
+ std::vector<double> val;
+ // val.resize(sol.size2()*sol.size1());
+ val.resize(sol.size1());
+ int counter = 0;
+ // for (int iC=0;iC<sol.size2();iC++)
+ printf("%g %g %g\n",sol(0,0),sol(1,0),sol(2,0));
+ for (int iR=0;iR<sol.size1();iR++)val[counter++] = sol(iR,0);
+ data[num] = val;
+ }
+ }
+
+ PView *pv = new PView (name, "ElementNodeData", _gm, data, 0.0, 1);
+ pv->getData()->writeMSH(name+".msh", false);
+ delete pv;
}
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index 879a6998b28c93e63064b221b6ebba6c681bdf52..567f37f5b90518b093623538b26bb4131da10b43 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -31,18 +31,19 @@ class dgSystemOfEquations {
int _order;
// dimension of the problem
int _dimension;
- // the solution
- fullMatrix<double> *_solution;
- // the right hand side (a temporary vector)
- fullMatrix<double> *_rightHandSide;
+ // solution and righ hand sides as a large arrays of doubles
+ int _dataSize;
+ double * _solution;
+ double * _rightHandSide;
+ // proxis of the solution and of the right hand side for each group;
+ std::vector<fullMatrix<double> *> _solutionProxys;
+ std::vector<fullMatrix<double> *> _rightHandSideProxys;
// groups of elements (volume terms)
std::vector<dgGroupOfElements*> _elementGroups;
// groups of faces (interface terms)
std::vector<dgGroupOfFaces*> _faceGroups;
// groups of faces (boundary conditions)
std::vector<dgGroupOfFaces*> _boundaryGroups;
- // sets up everything
- void setup(GModel *gm, dgConservationLaw *claw, int order);
dgSystemOfEquations(const dgSystemOfEquations &) {}
public:
// lua stuff
@@ -58,11 +59,15 @@ public:
int setConservationLaw (lua_State *L); // set the conservationLaw
int addBoundaryCondition (lua_State *L); // add a boundary condition : "physical name", "type", [options]
int setup (lua_State *L); // setup the groups and allocate
+ int L2Projection (lua_State *L); // assign the solution to a given function
dgSystemOfEquations(lua_State *L);
#endif // HAVE LUA
- dgSystemOfEquations(GModel *gm, dgConservationLaw *claw, int order);
+ inline const fullMatrix<double> getSolutionProxy (int iGroup, int iElement){
+ return fullMatrix<double> ( *_solutionProxys [iGroup] ,
+ iElement * _claw->nbFields(),_claw->nbFields());
+ }
+ void export_solution_as_is (const std::string &fileName);
~dgSystemOfEquations();
-
private:
};
diff --git a/Solver/luaFunction.cpp b/Solver/luaFunction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..74cde201e6848db93c4dde431028165f71154d1c
--- /dev/null
+++ b/Solver/luaFunction.cpp
@@ -0,0 +1,100 @@
+#include "luaFunction.h"
+#if defined(HAVE_LUA)
+#include <string>
+#include <vector>
+#include "function.h"
+// function that is defined in Lua
+
+class functionLua : public function {
+ lua_State *_L;
+ std::string _luaFunctionName;
+ std::vector<std::string> _dependenciesName;
+ int _nbCol;
+ private:
+ class data : public dataCacheDouble{
+ private:
+ dataCacheDouble &_uvw;
+ const functionLua *_function;
+ std::vector<dataCacheDouble *> _dependencies;
+ public:
+ data(const functionLua *f, dataCacheMap *m)
+ : _function(f),_uvw(m->get("UVW",this))
+ {
+ _dependencies.resize ( _function->_dependenciesName.size());
+ for (int i=0;i<_function->_dependenciesName.size();i++)
+ _dependencies[i] = &m->get(_function->_dependenciesName[i],this);
+ }
+ void _eval()
+ {
+ if(_value.size1()!=_uvw().size1() || _value.size2() != _function->_nbCol){
+ _value=fullMatrix<double>(_uvw().size1(),_function->_nbCol);
+ }
+ lua_getfield(_function->_L, LUA_GLOBALSINDEX, _function->_luaFunctionName.c_str());
+ for (int i=0;i< _dependencies.size();i++){
+ const fullMatrix<double> *data = &(*_dependencies[i])();
+ lua_pushlightuserdata (_function->_L, (void*) data);
+ }
+ lua_pushlightuserdata (_function->_L, &_value);
+ lua_call(_function->_L,_dependencies.size()+1,0); /* call Lua function */
+ }
+ };
+ public:
+ functionLua (int nbCol,
+ std::string &luaFunctionName,
+ std::vector<std::string> &dependenciesName,
+ lua_State *L)
+ : _luaFunctionName(luaFunctionName), _dependenciesName(dependenciesName),_L(L),_nbCol(nbCol)
+ {
+ }
+
+ dataCacheDouble *newDataCache(dataCacheMap *m)
+ {
+ return new data(this,m);
+ }
+};
+
+static int newLuaFunction (lua_State *L){
+ int n = lua_gettop(L);
+ std::vector<std::string> dependenciesName;
+ int nbFields = luaL_checkinteger(L, 1);
+ std::string luaFunctionName = std::string(luaL_checkstring(L, 2));
+ for (int i=3;i<=n;i++)
+ dependenciesName.push_back(std::string(luaL_checkstring(L, i)));
+ int ITER = 1;
+ std::string functionName = luaFunctionName;
+ while (function::get(functionName, true)){
+ char toto[256];
+ sprintf(toto,"%d",ITER);
+ functionName = luaFunctionName + toto;
+ }
+ function::add(functionName,new functionLua(nbFields,luaFunctionName,dependenciesName,L));
+ lua_pushstring(L, functionName.c_str());
+
+ return 1;
+}
+
+static int newConstantFunction (lua_State *L){
+ fullMatrix<double> * _ud = (fullMatrix<double> *) lua_touserdata(L, 1);
+
+ int ITER = 1;
+ std::string functionName = "ConstantFunction";
+ while (function::get(functionName, true)){
+ char toto[256];
+ sprintf(toto,"%d",ITER);
+ functionName = functionName + toto;
+ }
+
+ function::add(functionName,function::newFunctionConstant(*_ud));
+ lua_pushstring(L, functionName.c_str());
+ return 1;
+}
+
+int RegisterFunctions (lua_State *L) {
+ luaL_reg fcts[] = {{"lua", newLuaFunction },
+ {"constant", newConstantFunction }};
+ luaL_register(L, "createFunction", fcts);
+};
+
+
+
+#endif // HAVE LUA
diff --git a/Solver/luaFunction.h b/Solver/luaFunction.h
new file mode 100644
index 0000000000000000000000000000000000000000..bc161db373b85a9b23399b9a0b661c93ef3644e8
--- /dev/null
+++ b/Solver/luaFunction.h
@@ -0,0 +1,13 @@
+#include "GmshConfig.h"
+#ifndef _LUA_FUNCTION_H_
+#define _LUA_FUNCTION_H_
+#if defined(HAVE_LUA)
+// include lua stuff
+extern "C" {
+#include "lua.h"
+#include "lauxlib.h"
+#include "lualib.h"
+}
+int RegisterFunctions (lua_State *L);
+#endif // HAVE LUA
+#endif // _LUA_FUNCTION_H_
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index f0006bbbd46a268171e361ff834aacf25926d872..97becf7def0719b999d0b9f816d8108b24788fbd 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -47,7 +47,7 @@ class DI_Point
// assignment
DI_Point & operator=(const DI_Point & rhs);
// destructor
- virtual ~DI_Point () {Ls.clear();}
+ virtual ~DI_Point () {}
// add a levelset value (adjusted to 0 if ls<ZERO_LS_TOL) into the vector Ls
void addLs (const double ls);
// add a levelset value evaluated into the element e
@@ -218,6 +218,7 @@ class DI_Element
// domain elements : -1 = outside / +1 = inside
// interface elements : tag of the levelset that created the element
// -1 = out of the domain border
+ // DI_Point *pts_; // vertices
DI_Point **pts_; // vertices
DI_Point **mid_; // middle vertices
int polOrder_; // polynomial order of the shape functions
@@ -624,8 +625,10 @@ class DI_Quad : public DI_Element
class DI_Tetra : public DI_Element
{
+ //DI_Point pts_[4]; // vertices
public:
DI_Tetra () {
+ //pts_ = new DI_Point [4];
pts_ = new DI_Point*[4];
for(int i = 0; i < 4; i++) pts_[i] = NULL;
}
@@ -633,6 +636,7 @@ class DI_Tetra : public DI_Element
double x2, double y2, double z2, double x3, double y3, double z3)
{
pts_ = new DI_Point*[4];
+ //pts_[0] = DI_Point(x0, y0, z0);
pts_[0] = new DI_Point(x0, y0, z0);
pts_[1] = new DI_Point(x1, y1, z1);
pts_[2] = new DI_Point(x2, y2, z2);