diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index 0312dd4b40440f8fa33c684cb537c248deda1126..6b601c9d8a313430edf2ce3e2ecf73e91ccc9592 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -36,6 +36,18 @@ extern "C" {
std::complex<double> *alpha, std::complex<double> *a, int *lda,
std::complex<double> *x, int *incx, std::complex<double> *beta,
std::complex<double> *y, int *incy);
+ void F77NAME(daxpy)(int *n, double *alpha, double *x, int *incx, double *y, int *incy );
+
+
+}
+
+template<>
+void fullVector<double>::axpy(fullVector<double> &x,double alpha)
+{
+ // for (int i=0;i<_r;i++) _data[i] += alpha * x._data[i];
+ // return;
+ int M = _r, INCX = 1, INCY = 1;
+ F77NAME(daxpy)(&M, &alpha, x._data,&INCX, _data, &INCY);
}
template<>
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 4fe957c33d55f753f0c4718642a56924679b3d8c..da87ddb23f88060bde412b2ac51ecf56f4ee2a81 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -71,6 +71,14 @@ class fullVector
for(int i = 0; i < _r; ++i) s += _data[i] * other._data[i];
return s;
}
+ // y <- y + alpha * x
+ void axpy(fullVector<scalar> &x, scalar alpha=1.)
+#if !defined(HAVE_BLAS)
+ {
+ for (int i=0;i<_r;i++) _data[i] += alpha * x._data[i];
+ }
+#endif
+ ;
void print(const char *name="") const
{
printf("Printing vector %s:\n", name);
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 0203bcf63ad39d616755ab5258b3600845fdc2fd..da3d0ebfa66d2dfd1b5e66fb84f30aec8cc7aa2c 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -1,6 +1,7 @@
#include <set>
#include "dgAlgorithm.h"
#include "dgGroupOfElements.h"
+#include "dgSystemOfEquations.h"
#include "dgConservationLaw.h"
#include "GEntity.h"
#include "MElement.h"
@@ -217,9 +218,9 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
}
void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
- const dgGroupOfElements & group,
- fullMatrix<double> &residu,
- fullMatrix<double> &sol)
+ const dgGroupOfElements & group,
+ fullMatrix<double> &residu,
+ fullMatrix<double> &sol)
{
fullMatrix<double> residuEl;
fullMatrix<double> solEl;
@@ -234,59 +235,72 @@ 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
- std::vector<fullMatrix<double> *> &residu,
- std::vector<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
+ dgDofContainer &sol,
+ dgDofContainer &resd,
+ 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]);
+ // 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 nbFields = claw.nbFields();
+
+ dgDofContainer K (eGroups,claw);
+ dgDofContainer Unp (eGroups,claw);
+
+ K._data->scale(0.0);
+ K._data->axpy(*(sol._data));
+ Unp._data->scale(0.0);
+ Unp._data->axpy(*(sol._data));
+
+ for(int j=0; j<orderRK;j++){
+ if(j){
+ K._data->scale(b[j]);
+ K._data->axpy(*(sol._data));
}
- sol=Unp;
+ // printf("iter %d sol size = %d\n",j,sol._dataSize);
+ this->residual(claw,eGroups,fGroups,bGroups,K._dataProxys,resd._dataProxys);
+ K._data->scale(0.);
+ for(int k=0;k < eGroups.size();k++) {
+ int nbNodes = eGroups[k]->getNbNodes();
+ for(int i=0;i<eGroups[k]->getNbElements();i++) {
+ residuEl.setAsProxy(*(resd._dataProxys[k]),i*nbFields,nbFields);
+ KEl.setAsProxy(*(K._dataProxys[k]),i*nbFields,nbFields);
+ iMassEl.setAsProxy(eGroups[k]->getInverseMassMatrix(),i*nbNodes,nbNodes);
+ iMassEl.mult(residuEl,KEl);
+ }
}
-
- */
+ Unp._data->axpy(*(K._data),a[j]);
+ }
+
+ for (int i=0;i<sol._dataSize;i++){
+ // printf("tempSol[%d] = %g\n",i,(*tempSol._data)(i));
+ // memcp
+ (*sol._data)(i)=(*Unp._data)(i);
+ }
}
void dgAlgorithm::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
- )
+ const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &solutionLeft,
+ fullMatrix<double> &residual // residual !! at faces nodes
+ )
{
int nbFields = claw.nbFields();
const dgBoundaryCondition *boundaryCondition = claw.getBoundaryCondition(group.getBoundaryTag());
@@ -409,8 +423,8 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
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;
+ if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
+ if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
}
fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
@@ -424,8 +438,8 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
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;
+ if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
+ if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
}
fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index dfad3132e488ba949dd3a9674a8e4c4e08e69ca1..566cf906d0c268872dcebe91a3867f6284790bf9 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -7,6 +7,7 @@ class GModel;
class dgGroupOfElements;
class dgGroupOfFaces;
class dgConservationLaw;
+class dgDofContainer;
class dgTerm;
class dgAlgorithm {
@@ -44,8 +45,8 @@ class dgAlgorithm {
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,
+ dgDofContainer &residu,
+ dgDofContainer &sol,
int orderRK=4);
static void multAddInverseMassMatrix ( /*dofManager &dof,*/
const dgGroupOfElements & group,
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 6cc5af311124f95cf17de611f4f7a66a5cb24db8..acaa38827677a29cc3f6f91d85248e3f53f803f3 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -13,9 +13,7 @@ const char dgSystemOfEquations::className[] = "dgSystemOfEquations";
// Define the methods we will expose to Lua
#define _method(class, name) {#name, &class::name}
Luna<dgSystemOfEquations>::RegType dgSystemOfEquations::methods[] = {
- _method(dgSystemOfEquations, getRHS),
- _method(dgSystemOfEquations, computeRHS),
- _method(dgSystemOfEquations, getSolution),
+ _method(dgSystemOfEquations, RK44),
_method(dgSystemOfEquations, exportSolution),
_method(dgSystemOfEquations, setConservationLaw),
_method(dgSystemOfEquations, setOrder),
@@ -56,6 +54,7 @@ dgSystemOfEquations::dgSystemOfEquations(lua_State *L){
if (!obj)throw;
_gm = obj->getGModel();
_dimension = _gm->getNumRegions() ? 3 : 2;
+ _solution = 0;
}
// set the conservation law as a string (for now)
@@ -98,8 +97,8 @@ int dgSystemOfEquations::addBoundaryCondition (lua_State *L){
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]);
+ _algo->residualVolume(Law,*_elementGroups[i],*_solution->_dataProxys[i],*_rightHandSide->_dataProxys[i]);
+ _algo->multAddInverseMassMatrix(*_elementGroups[i],*_rightHandSide->_dataProxys[i],*_solution->_dataProxys[i]);
}
}
@@ -113,60 +112,35 @@ int dgSystemOfEquations::setup(lua_State *L){
_faceGroups,
_boundaryGroups);
// compute the total size of the soution
- _dataSize = 0;
- int nbFields = _claw->nbFields();
- 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;
- }
+ _solution = new dgDofContainer(_elementGroups,*_claw);
+ _rightHandSide = new dgDofContainer(_elementGroups,*_claw);
+ // printf("aaaaaaaaaaaaa\n");
}
-int dgSystemOfEquations::getSolution(lua_State *L){
- lua_pushlightuserdata (L, _solution);
- return 1;
-}
-
-int dgSystemOfEquations::computeRHS(lua_State *L){
- _algo->residual(*_claw, _elementGroups, _faceGroups, _boundaryGroups, _solutionProxys, _rightHandSideProxys);
- lua_pushlightuserdata (L, _solution);
- return 1;
-}
-
-int dgSystemOfEquations::getRHS(lua_State *L){
- lua_pushlightuserdata (L, _rightHandSide);
+int dgSystemOfEquations::RK44(lua_State *L){
+ double dt = luaL_checknumber(L, 1);
+ _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide);
+ double normSolution = _solution->_data->norm();
+ lua_pushnumber (L, normSolution);
return 1;
}
int dgSystemOfEquations::exportSolution(lua_State *L){
std::string outputFile(luaL_checkstring(L, 1));
export_solution_as_is(outputFile);
+ return 0;
}
#endif // HAVE_LUA
dgSystemOfEquations::~dgSystemOfEquations(){
for (int i=0;i<_elementGroups.size();i++){
- delete _solutionProxys[i];
- delete _rightHandSideProxys[i];
delete _elementGroups[i];
}
- delete [] _solution;
- delete [] _rightHandSide;
+ if (_solution){
+ delete _solution;
+ delete _rightHandSide;
+ }
}
void dgSystemOfEquations::export_solution_as_is (const std::string &name){
@@ -229,3 +203,31 @@ void dgSystemOfEquations::export_solution_as_is (const std::string &name){
pv->getData()->writeMSH(name+".msh", false);
delete pv;
}
+
+dgDofContainer::dgDofContainer (std::vector<dgGroupOfElements*> &elementGroups, const dgConservationLaw &claw){
+ _dataSize = 0;
+ int nbFields = claw.nbFields();
+ 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
+ _data = new fullVector<double>(_dataSize);
+ // create proxys for each group
+ int offset = 0;
+ _dataProxys.resize(elementGroups.size());
+ for (int i=0;i<elementGroups.size();i++){
+ int nbNodes = elementGroups[i]->getNbNodes();
+ int nbElements = elementGroups[i]->getNbElements();
+ _dataProxys[i] = new fullMatrix<double> (&(*_data)(offset*sizeof(double)),nbNodes, nbFields*nbElements);
+ offset += nbNodes*nbFields*nbElements;
+ }
+ // printf("datasize = %d\n",_dataSize);
+}
+
+dgDofContainer::~dgDofContainer (){
+ if (!_dataSize)return;
+ for (int i=0;i<_dataProxys.size();++i) delete _dataProxys[i];
+ delete _data;
+}
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index 567f37f5b90518b093623538b26bb4131da10b43..45c5108a6a9fe5cc363dba3b7ec0094dbabd993f 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -19,6 +19,18 @@ extern "C" {
#include "luna.h"
#endif // HAVE LUA
+struct dgDofContainer {
+private:
+ dgDofContainer (const dgDofContainer&);
+public:
+ int _dataSize; // the full data size i.e. concerning all groups
+ fullVector<double> * _data; // the full data itself
+ std::vector<fullMatrix<double> *> _dataProxys; // proxys
+ dgDofContainer (std::vector<dgGroupOfElements*> &groups, const dgConservationLaw &claw);
+ ~dgDofContainer ();
+};
+
+
class dgSystemOfEquations {
// the mesh and the model
GModel *_gm;
@@ -31,13 +43,9 @@ class dgSystemOfEquations {
int _order;
// dimension of the problem
int _dimension;
- // 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;
+ // solution and righ hand sides
+ dgDofContainer *_solution;
+ dgDofContainer *_rightHandSide;
// groups of elements (volume terms)
std::vector<dgGroupOfElements*> _elementGroups;
// groups of faces (interface terms)
@@ -51,19 +59,17 @@ public:
static const char className[];
static Luna<dgSystemOfEquations>::RegType methods[];
static void Register(lua_State *L);
- int computeRHS (lua_State *L); // get the Right Hand
- int getRHS (lua_State *L); // get the Right Hand
- int getSolution (lua_State *L); // get the Solution
int exportSolution (lua_State *L); // export the solution
int setOrder (lua_State *L); // set the polynomial order
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
+ int RK44 (lua_State *L); // assign the solution to a given function
dgSystemOfEquations(lua_State *L);
#endif // HAVE LUA
inline const fullMatrix<double> getSolutionProxy (int iGroup, int iElement){
- return fullMatrix<double> ( *_solutionProxys [iGroup] ,
+ return fullMatrix<double> ( *_solution->_dataProxys [iGroup] ,
iElement * _claw->nbFields(),_claw->nbFields());
}
void export_solution_as_is (const std::string &fileName);
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index 1d24a0a5c162c32c23bf8491f233d8733efbea87..2665373bb560ec30371424fc0182a9277a690a13 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -43,11 +43,11 @@ class DI_Point
DI_Point (double x, double y, double z, const double ls) : x_(x), y_(y), z_(z) {addLs(ls);}
DI_Point (double x, double y, double z, const DI_Element *e,
const std::vector<const gLevelset *> &RPNi) : x_(x), y_(y), z_(z) {computeLs(e, RPNi);}
- DI_Point(const DI_Point &p) : x_(p.x()), y_(p.y()), z_(p.z()) {Ls.clear(); Ls = p.Ls;}
+ // DI_Point(const DI_Point &p) : x_(p.x()), y_(p.y()), z_(p.z()) {Ls.clear(); Ls = p.Ls;}
// assignment
DI_Point & operator=(const DI_Point & rhs);
// destructor
- virtual ~DI_Point () {}
+ // 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