dg : organise files into modules

git-svn-id: https://geuz.org/svn/gmsh/trunk@16963 59c8cabf-10d6-4a28-b43f-3936b85956c4

CMakeLists.txt

@@ -119,115 +119,116 @@ if(ENABLE_LIM)
   endif (LIM_LIB)
 endif(ENABLE_LIM)
 
-include_directories(. conservationLaws slim ../../Common ../../Numeric ../../Geo ../../Mesh 
-  ../../Solver ../../Post ../../Plugin ../../Graphics 
-  ${GMSH_EXTERNAL_INCLUDE_DIRS} ${CMAKE_CURRENT_BINARY_DIR}/gmsh/Common 
-  ${CMAKE_CURRENT_BINARY_DIR} ${EXTERNAL_INCLUDES})
+include_directories(. conservationLaws slim ../../Common ../../Numeric ../../Geo ../../Mesh
+  ../../Solver ../../Post ../../Plugin ../../Graphics
+  ${GMSH_EXTERNAL_INCLUDE_DIRS} ${CMAKE_CURRENT_BINARY_DIR}/gmsh/Common
+  ${CMAKE_CURRENT_BINARY_DIR} ${EXTERNAL_INCLUDES}
+  conservationLaw dofContainer function limiter linearSystem mesh modules particle timeIntegrator tools)
 
 set(SRC
-  conservationLaws/dgConservationLawAdvection.cpp
-  conservationLaws/dgConservationLawAdvection2.cpp
-  conservationLaws/dgConservationLawAdvectionSUPG.cpp	
-  conservationLaws/dgConservationLawDiffusionSIP.cpp	
-  conservationLaws/dgConservationLawEulerAtm.cpp
-  conservationLaws/dgConservationLawFunction.cpp
-  conservationLaws/dgConservationLawLEE.cpp
-  conservationLaws/dgConservationLawMaxwell.cpp	
-  conservationLaws/dgConservationLawMaxwellExplicit.cpp
-  conservationLaws/dgConservationLawNavierStokes1d.cpp
-  conservationLaws/dgConservationLawNavierStokesIncomp2d.cpp	
-  conservationLaws/dgConservationLawNavierStokesIncomp3d.cpp
-  conservationLaws/dgConservationLawODE.cpp 
-  conservationLaws/dgConservationLawPerfectGas.cpp
-  conservationLaws/dgConservationLawPerfectGas2.cpp
-  conservationLaws/dgConservationLawPerfectGasALE.cpp
-  conservationLaws/dgConservationLawRichards.cpp
-  conservationLaws/dgConservationLawSW2DBench.cpp
-  conservationLaws/dgConservationLawSW3dContinuity.cpp
-  conservationLaws/dgConservationLawSW3dMomentum.cpp
-  conservationLaws/dgConservationLawSW3dTracer.cpp
-  conservationLaws/dgConservationLawSWnonInertia.cpp
-  conservationLaws/dgConservationLawShallowWater1d.cpp
-  conservationLaws/dgConservationLawShallowWater2d.cpp
-  conservationLaws/dgConservationLawShallowWaterAbsLayer2d.cpp
-  conservationLaws/dgConservationLawShallowWaterTracer1d.cpp
-  conservationLaws/dgConservationLawShallowWaterTracer2d.cpp
-  conservationLaws/dgConservationLawStokes3d.cpp
-  conservationLaws/dgConservationLawVerticalCoord.cpp
-  conservationLaws/dgConservationLawWave.cpp
-  conservationLaws/dgConservationLawWaveExplicit.cpp
-  conservationLaws/dgTerm.cpp
-  conservationLaws/dgTermControlVolumeFE.cpp
-  dgBlockMatrix.cpp
-  dgCGStructure.cpp
-  dgCircleTransform.cpp
-  dgConservationLaw.cpp
-  dgCoupling.cpp
-  dgCubedSphere.cpp
-  dgDofContainer.cpp
-  dgDofManager.cpp
-  dgElementVector.cpp
-  dgExtrusion.cpp
-  dgFilterBoydVandeven.cpp
-  dgFullMatrix.cpp
-  dgFunctionConstantByElement.cpp
-  dgFunctionEvaluator.cpp
-  dgFunctionGlobalOp.cpp
-  dgFunctionIntegrator.cpp
-  dgFunctionIntegratorInterface.cpp
-  dgFunctionLocalOp.cpp
-  dgGhostStructure.cpp
-  dgGroupOfElements.cpp
-  dgHOMesh.cpp
-  dgIntegrationMatrices.cpp
-  dgInterfaceVector.cpp
-  dgLimiter.cpp
-  dgLinearSystemBJacobi.cpp
-  dgLinearSystemExtrusion.cpp
-  dgLinearSystemMassMatrix.cpp
-  dgLinearSystemPETScShellDemo.cpp
-  dgMPI.cpp
-  dgMassAssembler.cpp
-  dgMesh.cpp
-  dgMesh2MeshProjection.cpp
-  dgMeshJacobian.cpp
-  dgMultirateERK.cpp
-  dgNeighbourMap.cpp
-  dgParticleTracker2D.cpp
-  dgSolver.cpp
-  dgSolverSupra.cpp
-  dgSpaceTransform.cpp
-  dgTagToFunction.cpp
-  dgTimer.cpp
-  dgTransformNodalValue.cpp
-  function.cpp
-  functionDerivator.cpp
-  functionGeneric.cpp
-  functionLookupTable.cpp
-  slim/dgEddyTransportFlux.cpp
-  slim/dgFEGradient.cpp
-  slim/dgGOTMWrapper.cpp 
-  slim/dgSW3dVerticalModel.cpp
-  slim/ncDataDz.cpp
-  slim/slim3dEquations.cpp
-  slim/slim3dExporter.cpp
-  slim/slim3dFunctions.cpp
-  slim/slim3dSolver.cpp
-  slim/slim3dTimeIntegrator.cpp
-  slim/slim3dUtils.cpp
-  slim/slimDate.cpp
-  slim/slimFields.cpp
-  slim/slimFunction.cpp
-  slim/slimGebco.cpp
-  slim/slimIceDyn.cpp
-  slim/slimIceDynSteady.cpp
-  slim/slimIterateHelper.cpp
-  slim/slimLonLat.cpp
-  slim/slimLim.cpp
-  slim/slimMovingBathWettingDrying.cpp
-  slim/slimStructData.cpp
-  slim/slimTemporalSerie.cpp
-  slim/slimTpxo.cpp
+  conservationLaw/dgConservationLaw.cpp
+  conservationLaw/dgConservationLawFunction.cpp
+  conservationLaw/dgCoupling.cpp
+  conservationLaw/dgMassAssembler.cpp
+  conservationLaw/dgTerm.cpp
+  dofContainer/dgDofContainer.cpp
+  dofContainer/dgGhostStructure.cpp
+  function/dgFunctionConstantByElement.cpp
+  function/dgFunctionEvaluator.cpp
+  function/dgFunctionGlobalOp.cpp
+  function/dgFunctionIntegrator.cpp
+  function/dgFunctionIntegratorInterface.cpp
+  function/dgFunctionLocalOp.cpp
+  function/dgTagToFunction.cpp
+  function/function.cpp
+  function/functionDerivator.cpp
+  function/functionGeneric.cpp
+  function/functionLookupTable.cpp
+  limiter/dgFilterBoydVandeven.cpp
+  limiter/dgLimiter.cpp
+  linearSystem/dgDofManager.cpp
+  linearSystem/dgLinearSystemBJacobi.cpp
+  linearSystem/dgLinearSystemExtrusion.cpp
+  linearSystem/dgLinearSystemMassMatrix.cpp
+  linearSystem/dgLinearSystemPETScShellDemo.cpp
+  mesh/dgCGStructure.cpp
+  mesh/dgCircleTransform.cpp
+  mesh/dgElementVector.cpp
+  mesh/dgExtrusion.cpp
+  mesh/dgGroupOfElements.cpp
+  mesh/dgHOMesh.cpp
+  mesh/dgInterfaceVector.cpp
+  mesh/dgMesh.cpp
+  mesh/dgMeshJacobian.cpp
+  mesh/dgNeighbourMap.cpp
+  mesh/dgSpaceTransform.cpp
+  modules/atmosphere/dgConservationLawEulerAtm.cpp
+  modules/atmosphere/dgCubedSphere.cpp
+  modules/generic/dgConservationLawAdvection.cpp
+  modules/generic/dgConservationLawAdvection2.cpp
+  modules/generic/dgConservationLawAdvectionSUPG.cpp
+  modules/generic/dgConservationLawDiffusionSIP.cpp
+  modules/generic/dgConservationLawODE.cpp
+  modules/lee/dgConservationLawLEE.cpp
+  modules/maxwell/dgConservationLawMaxwell.cpp
+  modules/maxwell/dgConservationLawMaxwellExplicit.cpp
+  modules/navierStokesIncomp/dgConservationLawNavierStokes1d.cpp
+  modules/navierStokesIncomp/dgConservationLawNavierStokesIncomp2d.cpp
+  modules/navierStokesIncomp/dgConservationLawNavierStokesIncomp3d.cpp
+  modules/navierStokesIncomp/dgConservationLawStokes3d.cpp
+  modules/perfectGas/dgConservationLawPerfectGas.cpp
+  modules/perfectGas/dgConservationLawPerfectGas2.cpp
+  modules/perfectGas/dgConservationLawPerfectGasALE.cpp
+  modules/shallowWater/dgConservationLawSW2DBench.cpp
+  modules/shallowWater/dgConservationLawShallowWater1d.cpp
+  modules/shallowWater/dgConservationLawShallowWater2d.cpp
+  modules/shallowWater/dgConservationLawShallowWaterAbsLayer2d.cpp
+  modules/shallowWater/dgConservationLawShallowWaterTracer1d.cpp
+  modules/shallowWater/dgConservationLawShallowWaterTracer2d.cpp
+  modules/shallowWater/slimMovingBathWettingDrying.cpp
+  modules/slim/dgConservationLawSW3dContinuity.cpp
+  modules/slim/dgConservationLawSW3dMomentum.cpp
+  modules/slim/dgConservationLawSW3dTracer.cpp
+  modules/slim/dgConservationLawVerticalCoord.cpp
+  modules/slim/dgEddyTransportFlux.cpp
+  modules/slim/dgFEGradient.cpp
+  modules/slim/dgGOTMWrapper.cpp
+  modules/slim/dgSW3dVerticalModel.cpp
+  modules/slim/ncDataDz.cpp
+  modules/slim/slim3dEquations.cpp
+  modules/slim/slim3dExporter.cpp
+  modules/slim/slim3dFunctions.cpp
+  modules/slim/slim3dSolver.cpp
+  modules/slim/slim3dTimeIntegrator.cpp
+  modules/slim/slim3dUtils.cpp
+  modules/slim/slimDate.cpp
+  modules/slim/slimFields.cpp
+  modules/slim/slimFunction.cpp
+  modules/slim/slimGebco.cpp
+  modules/slim/slimIceDyn.cpp
+  modules/slim/slimIceDynSteady.cpp
+  modules/slim/slimIterateHelper.cpp
+  modules/slim/slimLim.cpp
+  modules/slim/slimLonLat.cpp
+  modules/slim/slimStructData.cpp
+  modules/slim/slimTemporalSerie.cpp
+  modules/slim/slimTpxo.cpp
+  modules/subsurface/dgConservationLawRichards.cpp
+  modules/subsurface/dgConservationLawSWnonInertia.cpp
+  modules/subsurface/dgTermControlVolumeFE.cpp
+  modules/supra/dgSolverSupra.cpp
+  modules/wave/dgConservationLawWave.cpp
+  modules/wave/dgConservationLawWaveExplicit.cpp
+  particle/dgParticleTracker2D.cpp
+  timeIntegrator/dgMultirateERK.cpp
+  timeIntegrator/dgSolver.cpp
+  timeIntegrator/dgTransformNodalValue.cpp
+  tools/dgBlockMatrix.cpp
+  tools/dgFullMatrix.cpp
+  tools/dgIntegrationMatrices.cpp
+  tools/dgMPI.cpp
+  tools/dgMesh2MeshProjection.cpp
+  tools/dgTimer.cpp
 )
 set_property(TARGET gmsh APPEND PROPERTY EXCLUDE_FROM_ALL ON)
 

conservationLaw/dgConservationLaw.cpp

+#include "dgConservationLaw.h"
+#include "dgGroupOfElements.h"
+#include "dgDofContainer.h"
+
+void dgConservationLaw::computeAllTerms(double t, dgDofContainer &solution, dgDofContainer &residual, bool scatterSolutionDof) const
+{
+  if ( scatterSolutionDof )
+     solution.scatterBegin();
+  residual.setAll(0.);
+  const dgGroupCollection &groups = *residual.getGroups();
+  for (int i=0; i < groups.getNbElementGroups(); i++) {
+    const dgGroupOfElements &group = *groups.getElementGroup(i);
+    computeElementTerms(group, t, solution, residual);
+  }
+  for(int i=0;i < groups.getNbFaceGroups() ; i++) {
+    dgGroupOfFaces &group = *groups.getFaceGroup(i);
+    bool oneGhost = false;
+    for (int j = 0; j < group.nConnection(); ++ j) {
+      oneGhost |= (group.elementGroup(j).getGhostPartition() >= 0);
+    }
+    if (!oneGhost)
+      computeFaceTerms(group, t, solution, residual);
+  }
+  if ( scatterSolutionDof )
+     solution.scatterEnd();
+  for(int i = 0; i < groups.getNbFaceGroups(); i++) {
+    dgGroupOfFaces &group = *groups.getFaceGroup(i);
+    bool oneGhost = false;
+    for (int j = 0; j < group.nConnection(); ++ j) {
+      oneGhost |= (group.elementGroup(j).getGhostPartition() >= 0);
+    }
+    if (oneGhost){
+      computeFaceTerms(group, t, solution, residual);
+    }
+  }
+}
+

conservationLaw/dgConservationLaw.h

+#ifndef _DG_CONSERVATION_LAW_H_
+#define _DG_CONSERVATION_LAW_H_
+
+class dgDofContainer;
+class dgDofManager;
+class dgGroupOfElements;
+class dgGroupOfFaces;
+
+//todo dgHOMesh should maybe not be exposed outside dgDofManager
+class dgHOMesh;
+template<class t> class dgFullMatrix;
+
+#include <string>
+#include <vector>
+
+//TODO : function and dgExtrusion should not be needed
+class function;
+class dgExtrusion;
+class dgMassAssembler;
+class dgConservationLaw {
+  public :
+  enum imexMode {IMEX_IMPLICIT, IMEX_ALL};
+  enum massFactorType {MASS_FACTOR_NO, MASS_FACTOR_CONSTANT, MASS_FACTOR_TIME_DEPENDENT};
+  virtual int getNbFields() const = 0;
+  virtual bool isConstantJac() const {return false;}
+  virtual bool isLinear() const {return false;}
+  virtual void setTimeDirection(int timeDir) {}
+  virtual int isTemporal(int iField) const {return 1;}
+  virtual void setImexMode(imexMode mode) {}
+  virtual massFactorType haveMassFactor(const dgGroupOfElements&) const {return MASS_FACTOR_NO;}
+  // todo access to dgResidual should be done through dgDofManager
+  virtual void computeElementTerms(const dgGroupOfElements &group, double t, dgDofContainer &solution, dgDofContainer &residual, dgDofManager *jacDof = NULL, dgFullMatrix<double> *mass = NULL, dgDofContainer *variable = NULL) const = 0;
+  virtual void computeElementTermsWithMass(const dgGroupOfElements &group, double t, dgDofContainer &solution, dgDofContainer &residual, dgDofManager *jacDof = NULL, dgMassAssembler *massAssembler = NULL, dgDofContainer *variable = NULL) const = 0;
+  virtual void computeFaceTerms(const dgGroupOfFaces &faces, double t, dgDofContainer &solutionDof, dgDofContainer &residual,  dgDofManager *jacDof = NULL, dgDofContainer *variable = NULL) const = 0;
+  virtual void computeFaceTermsWithMass(const dgGroupOfFaces &faces, double t, dgDofContainer &solutionDof, dgDofContainer &residual,  dgDofManager *jacDof = NULL, dgMassAssembler *massAssembler = NULL, dgDofContainer *variable = NULL) const = 0;
+  virtual void computeElementaryTimeSteps (const dgGroupOfElements & group, dgDofContainer *solution, std::vector<double> & DT, dgExtrusion *extrusion = NULL) const = 0;
+  void computeAllTerms(double t, dgDofContainer &solution, dgDofContainer &residual, bool scatterSolutionDof = true) const;
+  virtual void fixStrongBC(dgDofManager &dof, dgDofContainer *oldSolution = NULL, double alpha = 1.) const {};
+  virtual ~dgConservationLaw() {};
+
+  //todo this should be here
+  //virtual void fillSparsityPattern(dgDofManager *dof, bool hasInterfaceTerm) const = 0;
+  //and we can remove this
+  virtual bool haveFaceTerm(const dgGroupOfFaces &faces) const {return true;}
+
+  //todo everything bellow still relies on functions (but is optional) => should be adapted:
+  //for dgLimiter
+  virtual const function *getOutsideValue(const dgGroupOfFaces*, const function *in) const {return NULL;}
+  virtual const function *getClipToPhysics(const std::string tag) const {return NULL;}
+};
+
+#endif

conservationLaw/dgConservationLawFunction.h

+#ifndef _DG_CONSERVATION_LAW_FUNCTION_H_
+#define _DG_CONSERVATION_LAW_FUNCTION_H_
+
+#include "dgConservationLaw.h"
+
+class dataCacheMap;
+class function;
+class dgConservationLawFunction;
+class dgDofContainer;
+class dgDofManager;
+class dgBoundaryCondition;
+class dgTerm;
+class dgTermAssembler;
+class dgFaceTermAssembler;
+class dgInterfacePsiTerm;
+class dgIntegrationMatrices;
+class dgExtrusion;
+class function;
+template <class t> class fullMatrix;
+
+
+#include "GmshMessage.h"
+#include <map>
+
+/**
+ * A conservation law is defined a convective flux (f), a diffusive flux(g), a source term(r)
+ * and a set of boundary conditions.\n
+ * \f[\partial_t L(u) = \nabla \cdot (\mathbf{f}(u,\nabla u,forcings)) 
+ * + r(u,forcings))\f]
+ */
+class dgStrongBoundaryCondition{
+  public:
+  int dim;
+  std::string physicalTag;
+  const function *f;
+  dgStrongBoundaryCondition (int dim_, const std::string physicalTag_, const function *f_): dim(dim_),physicalTag(physicalTag_), f(f_){};
+};
+
+class dgConservationLawFunction :public dgConservationLaw{
+ public:
+  typedef std::map<const std::string, const function*> termMap;
+ private:
+  std::vector<dgStrongBoundaryCondition> _strongBoundaryConditions;
+  typedef std::map<const std::string, dgBoundaryCondition*> bcMap;
+  bcMap _boundaryConditions;
+  bool _built; //cfr checkSetup
+ protected :
+  std::vector<int> temporal; //0: untemporal; 1:temporal; 2:mass-lumped
+  int _nbf;
+  termMap _volumeTerm0;    // DG-term for source
+  termMap _volumeTerm1;    // DG-term for convective and diffusive flux
+  termMap _interfaceTerm0;
+  termMap _interfaceTerm1N;
+  termMap _bifurcationTerm0;
+  termMap _bifurcationTerm1N;
+  termMap _maximumConvectiveSpeed;
+  termMap _diffusivity;
+  termMap _clipToPhysics;
+  std::map<const std::string, std::pair<const function*, bool> > _massFactor; // the boolean indicates if the term is constant in time
+  static inline const function * functionForTag(const termMap &map, const std::string tag)
+  {
+    termMap::const_iterator it = map.find(tag);
+    if (it == map.end()) it = map.find("");
+    if (it == map.end()) return NULL;
+    return it->second;
+  }
+ public:
+
+  int isTemporal(int iField) const { return temporal[iField]; }
+
+  // Constructor - Setup - Destructor
+  dgConservationLawFunction(int nbf) {
+    _nbf = nbf;
+    _built = false;
+    temporal.resize(nbf);
+    for (size_t k = 0; k < temporal.size(); k++) temporal[k] = true;
+  }
+  virtual void setup() {}
+  virtual ~dgConservationLawFunction() {}
+
+  // Check setup
+  inline void checkSetup () const {if(!_built) { const_cast<dgConservationLawFunction*>(this)->setup(); const_cast<dgConservationLawFunction*>(this)->_built = true;}}
+  inline void checkNonSetup () {if(_built) Msg::Fatal("this conservation law is already _built, you cannot set this parameter now");}
+
+  // Get information
+  /** Return the number of fields composing the unknowns of this conservation law.
+   * For vectorial fields, each components is counted as one field.
+   */
+  int getNbFields() const {return _nbf;}
+  const function *getClipToPhysics(const std::string tag) const;
+  const function *getInterfaceTerm0(const std::string tag) const;
+  const function *getdiffusivity(const std::string tag) const;
+  const function *getMaximumConvectiveSpeed(const std::string tag) const;
+
+  /** Associate a mesh tag with a boundary condition.*/
+  void addBoundaryCondition(std::string tag, dgBoundaryCondition *condition);
+  void addBoundaryCondition(const std::vector<std::string> tags, dgBoundaryCondition *condition);
+  dgBoundaryCondition *getBoundaryCondition(const std::string tag) const;
+  virtual const function *getOutsideValue(const dgGroupOfFaces*, const function *in) const;
+  bool haveFaceTerm(const dgGroupOfFaces &faces) const;
+
+  /** Add a strong BC to a physical entity.
+   * Give tag and dim and constant function (value_field0, value_field1,..., flag_field0, flag_field1,...). 
+   * Fields with flags >0 are fixed
+   * Strong boundary condition are not implemented in Discontinuous Galerkin
+   */
+  void addStrongBoundaryCondition (int dim, const std::string tag, const function *f);
+
+  /**
+   * Create a new boundary condition based on external data.
+   * The interface term with the provided tag computes the fluxes using the provided values as external solution");
+  */
+  dgBoundaryCondition *newOutsideValueBoundary(const std::string tag, const function *outsideValueFunction);
+  dgBoundaryCondition *newOutsideValueBoundary(const function *outsideValueFunction)
+    {return newOutsideValueBoundary("", outsideValueFunction);}
+
+  /**
+   * Create a new boundary condition based on external data. Generic version
+   * The interface term with the provided tag computes the fluxes using the provided values as external solution");
+  */
+  dgBoundaryCondition *newOutsideValueBoundaryGeneric(const std::string tag, const std::vector<const function*> toReplace,
+                                                                             const std::vector<const function*> replaceBy);
+  dgBoundaryCondition *newOutsideValueBoundaryGeneric(const std::vector<const function*> toReplace, const std::vector<const function*> replaceBy)
+    {return newOutsideValueBoundaryGeneric("", toReplace, replaceBy);}
+
+  /**
+   * Create a new boundary condition based on external data. Generic version, with the 2 which could depend on the 1
+   * The interface term with the provided tag computes the fluxes using the provided values as external solution");
+  */
+  dgBoundaryCondition *newOutsideValueBoundaryGeneric(const std::string tag,
+                                                      const std::vector<const function*> toReplace1,
+                                                      const std::vector<const function*> replaceBy1,
+                                                      const std::vector<const function*> toReplace2,
+                                                      const std::vector<const function*> replaceBy2);
+  dgBoundaryCondition *newOutsideValueBoundaryGeneric(const std::string tag,
+                                                      const std::map<const function*, const function*> *map,
+                                                      const std::map<const function*, const function*> *map2=NULL);
+  /**
+   * Create a new boundary condition with a given flux (no other fluxes will be computed)
+   */
+  dgBoundaryCondition *newNeumannBoundary(const function *fluxFunction);
+
+  /**
+  * Create a new boundary condition which set to 0 all the fluxes through the interfaces.
+  */
+  dgBoundaryCondition *new0FluxBoundary(const std::string tag="");
+
+  /**
+  * Create a new boundary condition using the values computed inside the domain as boundary values. The interface term with the provided tag computes the fluxes using the internal values as internal AND external values.
+  */
+  dgBoundaryCondition *newSymmetryBoundary(const std::string tag = "");
+
+public:
+  virtual dgTermAssembler *getTerm(const dgGroupOfElements &group) const;
+  virtual dgFaceTermAssembler *getFaceTerm(const dgGroupOfFaces &group) const;
+
+  static double muFactor(dataCacheMap *m);
+  void computeElementTerms(const dgGroupOfElements &group, double t, dgDofContainer &solution, dgDofContainer &residual, dgDofManager *jacDof = NULL, dgFullMatrix<double> *mass = NULL, dgDofContainer *variable = NULL) const;
+  virtual void computeElementTermsWithMass(const dgGroupOfElements &group, double t, dgDofContainer &solution, dgDofContainer &residual, dgDofManager *jacDof = NULL, dgMassAssembler *massAssembler = NULL, dgDofContainer *variable = NULL) const;
+  void computeFaceTerms(const dgGroupOfFaces &faces, double t, dgDofContainer &solutionDof, dgDofContainer &residual,  dgDofManager *jacDof = NULL, dgDofContainer *variable = NULL) const;
+  virtual void computeFaceTermsWithMass(const dgGroupOfFaces &faces, double t, dgDofContainer &solutionDof, dgDofContainer &residual,  dgDofManager *jacDof = NULL, dgMassAssembler *massAssembler = NULL, dgDofContainer *variable = NULL) const;
+  void computeElementaryTimeSteps (const dgGroupOfElements & group, // the group
+                                          dgDofContainer *solution, // the solution 
+                                          std::vector<double> & DT, // elementary time steps
+                                          dgExtrusion *extrusion = NULL) const;
+  /**Return the maximum time step in order to ensure stability */
+  double getMinOfTimeSteps(dgDofContainer *solution, dgExtrusion *extrusion = NULL) const;
+  void fixStrongBC(dgDofManager &dof, dgDofContainer *oldSolution = NULL, double alpha = 1.) const;
+  const function * getMassFactor(const std::string tag, bool &constMass) const;
+  massFactorType haveMassFactor(const dgGroupOfElements &group) const;
+};
+
+
+/**
+ * A weak boundary condition. 
+ * Boundary conditions should be associated with tag using dgConservationLawFunction::addBoundaryCondition.
+ */
+class dgBoundaryCondition {
+  protected :
+  dgConservationLawFunction *_claw;
+  const function *_term0, *_term1N;
+  std::string _type;
+  std::map<const function*,const function*> _outsideValue;
+ public: 
+
+  // Constructor - Destructor
+  dgBoundaryCondition () {
+    _term0 = NULL;
+    _term1N = NULL;
+    _type = "undefined";
+  }
+  virtual ~dgBoundaryCondition () {}
+
+  // Get information
+  virtual const function * getTerm0 () { return _term0; }
+  virtual const function * getTerm1N () { return _term1N; }
+  virtual const function * getOutsideValue (const function* f) { return _outsideValue[f]; } 
+  inline std::string getType() const { return _type; }
+};
+
+/*==============================================================================
+ * IP Term
+ *============================================================================*/
+
+function *dgNewIpTerm (int nbFields, const function *diffusiveFlux, const function *diffusivity, bool _useSip=true);
+function *dgNewIpTermStereo (int nbFields, const function *diffusiveFlux, const function * nuhh, const function * nuvv, double R, bool verticalOnly);
+function *dgNewIpTermAnisotropic(int nbFields, const function *diffusiveFlux, const function *ipPenalty);
+function *dgNewIpTermErnRef(int nbFields, const function *diffusiveFlux,const function *ipPenalty);
+function *dgNewSymmetricIpTerm (int nbFields, const function *diffusivity);
+function *dgNewSymIpTerm3d(int nbFields, const function *diffusivity);
+#endif
+
+
+
+
+
+
+
+
+
+

conservationLaw/dgCoupling.cpp

+#include "dgCoupling.h"
+#include "dgIntegrationMatrices.h"
+#include <algorithm>
+
+class dgCoupling_SW1D_SW2D::computeToIntegrate1:public function{
+  public:
+  fullMatrix<double> _solution;
+  computeToIntegrate1(const function *solution):function(3){
+    setArgument(_solution, solution);
+  }
+  void call(dataCacheMap *m, fullMatrix<double> &val){
+    for( int i=0; i<val.size1(); i++){
+      val(i,0) = _solution(i,0);
+      val(i,1) = _solution(i,1); 
+      val(i,2) = 0;
+    }
+  }
+};
+function *dgCoupling_SW1D_SW2D::computeToIntegrate1New(const function *solution){
+  return new computeToIntegrate1(solution);
+}
+
+class dgCoupling_SW1D_SW2D::computeToIntegrate2:public function{
+  public:
+  fullMatrix<double> _solution, _bath, normals;
+  computeToIntegrate2(const function *solution, const function *bath):function(3){
+    setArgument(_solution, solution);
+    setArgument(_bath, bath);
+    setArgument(normals,function::getNormals());
+  }
+  void call(dataCacheMap *m, fullMatrix<double> &val){
+    for( int i=0; i<val.size1(); i++){
+      val(i,0) = _solution(i, 0) + _bath(i, 0);
+      val(i,1) = (_solution(i, 1) * normals(i,0) + _solution(i, 2) * normals(i,1)) / val(i, 0); 
+      val(i,2) = (_solution(i, 1) * (-normals(i,1)) + _solution(i, 2) * normals(i,0)) / val(i, 0); 
+    }
+  }
+};
+function *dgCoupling_SW1D_SW2D::computeToIntegrate2New(const function *solution, const function *bath){
+  return new computeToIntegrate2(solution, bath);
+}
+
+class dgCoupling_SW1D_SW2D::computeOutsideSolution1:public function{
+  public:
+  fullMatrix<double> _integrals;
+  double _section0, _width;
+  computeOutsideSolution1(const function *integrals, const double section0, const double width):function(3){
+    setArgument(_integrals, integrals);
+    _section0 = section0;
+    _width = width;
+  }
+  void call(dataCacheMap *m, fullMatrix<double> &val){
+    for( int i=0; i<val.size1(); i++){
+      val(i,0) = (_integrals(0, 0) - _section0) / _width; // = elevation = (S - int h dx) /   
+      val(i,1) = _integrals(0, 1); //velocity (integral on one point change nothing)
+      val(i,2) = 0;