dgConservationLawSW3dContinuity.cpp 17.6 KB
Newer Older
1
#include "dgConfig.h"
lambrechts's avatar
lambrechts committed
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#include "dgConservationLawFunction.h"
#include "dgIntegrationMatrices.h"
#include "fullMatrix.h"
#include "dgGroupOfElements.h"
#include "dgDofContainer.h"
#include "function.h"
#include "dgConservationLawSW3dContinuity.h"
#include "Numeric.h"
#include <limits.h>
#include "slim3dFunctions.h"
#include "dgMesh.h"
#include "dgMeshJacobian.h"
#include "GmshDefines.h"
#include "dgTerm.h"
lambrechts's avatar
lambrechts committed
16
#include "functorMember.h"
17
#include "dgExtrusion.h"
lambrechts's avatar
lambrechts committed
18

lambrechts's avatar
lambrechts committed
19
20
21
void dgConservationLawSW3dContinuity::advection(functorCache &cache, fullMatrix<double> &val) const
{
  size_t nQP = cache.nEvaluationPoint();
22
  val.resize(nQP, 3, true);
lambrechts's avatar
lambrechts committed
23
24
  const fullMatrix<double> &uv = cache.get(*_uv);
  for(size_t i=0; i< nQP; i++) {
25
26
    val(i,0) = uv(i,0);
    val(i,1) = uv(i,1);
27
    //val(i,2) = 0;    //already set to zero in resize
lambrechts's avatar
lambrechts committed
28
  }
lambrechts's avatar
lambrechts committed
29
}
lambrechts's avatar
lambrechts committed
30

lambrechts's avatar
lambrechts committed
31
32
void dgConservationLawSW3dContinuity::riemann(functorCache &cache, fullMatrix<double> &val) const
{
33
  val.resize(cache.nEvaluationPoint(), 2, false);
lambrechts's avatar
lambrechts committed
34
35
36
37
38
  const fullMatrix<double> &normals = cache.get(*function::getNormals(), 0);
  const fullMatrix<double> &etaIn = cache.get(*_eta, 0);
  const fullMatrix<double> &etaOut = cache.get(*_eta, 1);
  const fullMatrix<double> &UIn = cache.get(*_uv, 0);
  const fullMatrix<double> &UOut = cache.get(*_uv, 1);
39
40
  const fullMatrix<double> &wOldIn  = cache.get(*_wOld, 0);
  const fullMatrix<double> &wOldOut = cache.get(*_wOld, 1);
41
  const fullMatrix<double> &bath = cache.get(*_bathymetryMinCG, 0);
lambrechts's avatar
lambrechts committed
42
43
44
45
46
47
48
49
50
51
  double g = slim3dParameters::g;
  // determine which element is below
  bool outIsTop = (normals(0, 2) > 0);
  for(int i=0; i< val.size1(); i++) {
    double nx = normals(i,0);
    double ny = normals(i,1);
    double nz = normals(i,2);
    double unIn = nx*UIn(i,0) + ny*UIn(i,1);
    double unOut = nx*UOut(i,0) + ny*UOut(i,1);
    double un = (unIn+unOut)/2;
52
53
54
55
56
57
58
59
60
61
    if (fabs(nz) < 1e-8) { // lateral face
      if ( _laxFriedrichsFactor > 1e-8 )
        un = ( unIn+unOut + _laxFriedrichsFactor* (wOldIn(i,0) - wOldOut(i,0)) )/2;
      else{
        double hMin = bath(i,0);
        double etaM = 0.5*(etaIn(i,0)+etaOut(i,0));
        double sq_g_h = sqrt(g/(hMin+etaM));
        un = (unIn+unOut + (etaIn(i,0)-etaOut(i,0))*sq_g_h)/2;
      }
    } else { // top-bottom face
lambrechts's avatar
lambrechts committed
62
63
64
65
      if (outIsTop == _propagateFromTop) {
        un = unOut;
      } else {
        un = unIn;
lambrechts's avatar
lambrechts committed
66
67
      }
    }
68
69
    val(i,0) = -un;
    val(i,1) = -val(i,0);
lambrechts's avatar
lambrechts committed
70
  }
lambrechts's avatar
lambrechts committed
71
}
lambrechts's avatar
lambrechts committed
72
73
74
75
76

dgConservationLawSW3dContinuity::dgConservationLawSW3dContinuity( int propagateFromTop ) :
  dgConservationLawFunction(1), _propagateFromTop(propagateFromTop){
    _uv = NULL;
    _eta = NULL;
77
    _wOld = NULL;
lambrechts's avatar
lambrechts committed
78
79
}
void dgConservationLawSW3dContinuity::setup() {
80
  _volumeTerm0[""] = NULL;
lambrechts's avatar
lambrechts committed
81
82
  _volumeTerm1[""] = newFunctorMember(*this, &dgConservationLawSW3dContinuity::advection);
  _interfaceTerm0[""] = newFunctorMember(*this, &dgConservationLawSW3dContinuity::riemann);
lambrechts's avatar
lambrechts committed
83
84
85
86
  if (!_uv)
    Msg::Fatal("dgConservationLawSW3dContinuity: uv function must be set");
  if (!_eta)
    Msg::Fatal("dgConservationLawSW3dContinuity: eta function must be set");
87
88
  if (!_wOld)
    Msg::Fatal("dgConservationLawSW3dContinuity: wOld function must be set");
lambrechts's avatar
lambrechts committed
89
90
91
92
93
94
}
dgConservationLawSW3dContinuity::~dgConservationLawSW3dContinuity() {
  if (_volumeTerm0[""]) delete _volumeTerm0[""];
  if (_volumeTerm1[""]) delete _volumeTerm1[""];
  if (_interfaceTerm0[""]) delete _interfaceTerm0[""];
}
lambrechts's avatar
lambrechts committed
95
96

void dgConservationLawSW3dContinuity::boundaryWall(functorCache &cache, fullMatrix<double> &val) const {
97
  val.resize(cache.nEvaluationPoint(), 1, false);
lambrechts's avatar
lambrechts committed
98
99
100
101
102
  const fullMatrix<double> &uv = cache.get(*_uv, 0);
  const fullMatrix<double> &normals = cache.get(*function::getNormals(), 0);
  cache.pushState();
  cache.setDefaultToSide0(true);
  fullMatrix<double> &uvExt = cache.set(*_uv, 1);
103
  uvExt.resize(cache.nEvaluationPoint(), 2, false);
lambrechts's avatar
lambrechts committed
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
  for(int i=0;i<val.size1(); i++) {
    double nx = normals(i,0);
    double ny = normals(i,1);
    double nz = normals(i,2);
    // compute normal/tangential velocity
    double un =  uv(i,0)*nx + uv(i,1)*ny;
    //double ut = -uv(i,0)*ny + uv(i,1)*nx;
    if (fabs(nz)>1e-8) { // non used!! (0 flux used for bottom)
      // symmetry for horizontal faces
      uvExt(i,0) = uv(i,0);
      uvExt(i,1) = uv(i,1);
    } else {
      // flip normal component
      uvExt(i, 0) = uv(i,0) - 2*un*nx;
      uvExt(i, 1) = uv(i,1) - 2*un*ny;
lambrechts's avatar
lambrechts committed
119
    }
lambrechts's avatar
lambrechts committed
120
121
122
    // symmetry for velocity and elevation
    //           externalUV(i,0) = _uv(i,0);
    //           externalUV(i,1) = _uv(i,1);
lambrechts's avatar
lambrechts committed
123
  }
lambrechts's avatar
lambrechts committed
124
125
126
127
128
129
  const fullMatrix<double> &lawTerm = cache.get(*getInterfaceTerm0(""));
  for(int i=0;i<val.size1(); i++) {
    val(i, 0) = lawTerm(i, 0);
  }
  cache.popState();
}
lambrechts's avatar
lambrechts committed
130
131

dgBoundaryCondition *dgConservationLawSW3dContinuity::newBoundaryWall(){
lambrechts's avatar
lambrechts committed
132
  return newBoundaryConditionMember(*this, &dgConservationLawSW3dContinuity::boundaryWall);
lambrechts's avatar
lambrechts committed
133
134
}

135
136
void dgConservationLawSW3dPressure::source(functorCache &cm, fullMatrix<double> &val) const {
    const fullMatrix<double> &rhoGrad = cm.get(*_rhoGrad); 
137
    double rho0 = slim3dParameters::rho0;
138
    int nQP = cm.nEvaluationPoint();
139
    val.resize(nQP,3, true);
140
    for(int i=0; i<nQP; i++) {
141
      // d /dx
142
      val(i,0) = rhoGrad(i,0)/rho0;
143
      // d /dy
144
      val(i,1) = rhoGrad(i,1)/rho0;
145
      // d /dz
146
      //val(i,2) = 0;     //already set to zero in resize
147
    }
148
149
150
151
}

void dgConservationLawSW3dPressure::advection(functorCache &cm, fullMatrix<double> &val) const {
    const fullMatrix<double> &rho = cm.get(*_rho);
lambrechts's avatar
lambrechts committed
152
    double rho0 = slim3dParameters::rho0;
153
    int nQP = cm.nEvaluationPoint();
154
    val.resize(nQP,6, true);
155
    for(int i=0; i<nQP; i++) {
lambrechts's avatar
lambrechts committed
156
      // d /dx
157
      val(i,0) = rho(i,0)/rho0;
158
      //val(i,1) = 0;    //already set to zero in resize
lambrechts's avatar
lambrechts committed
159
      // d /dy
160
      //val(i,2) = 0;
161
      val(i,3) = rho(i,0)/rho0;
lambrechts's avatar
lambrechts committed
162
      // d /dz
163
164
      //val(i,4) = 0;
      //val(i,5) = 0;
lambrechts's avatar
lambrechts committed
165
    }
166
167
168
169
170
171
172
}

void dgConservationLawSW3dPressure::riemann(functorCache &cm, fullMatrix<double> &val) const {
    const fullMatrix<double> &normals = cm.get(*function::getNormals(),0);
    const fullMatrix<double> &rhoIn = cm.get(*_rho,0);
    const fullMatrix<double> &rhoOut = cm.get(*_rho,1);
    double rho0 = slim3dParameters::rho0; 
lambrechts's avatar
lambrechts committed
173
    // determine which element is below
174
    bool outIsTop = normals(0, 2) > 0;
175
    int nQP = cm.nEvaluationPoint();
176
    val.resize(nQP, 2, false);
177
    for(int i=0; i<nQP; i++) {
lambrechts's avatar
lambrechts committed
178
179
180
      double nx = normals(i,0);
      double ny = normals(i,1);
      double nz = normals(i,2);
181
      double rho = 0;
lambrechts's avatar
lambrechts committed
182
      if (fabs(nz) < 1e-6) { // horizontal -> vertical face
183
        rho = 0.5*(rhoIn(i,0)+rhoOut(i,0));
lambrechts's avatar
lambrechts committed
184
185
      } else { // horizontal face, top/bottom
        if (outIsTop == _propagateFromTop) {
186
          rho = rhoOut(i,0);
lambrechts's avatar
lambrechts committed
187
        } else {
188
          rho = rhoIn(i,0);
lambrechts's avatar
lambrechts committed
189
190
        }
      }
191
192
      val(i,0) = - (rho/rho0*nx);
      val(i,1) = - (rho/rho0*ny);
lambrechts's avatar
lambrechts committed
193
    }
194
}
lambrechts's avatar
lambrechts committed
195

196
dgConservationLawSW3dPressure::dgConservationLawSW3dPressure(int propagateFromTop): dgConservationLawFunction(2), _propagateFromTop(propagateFromTop) {
lambrechts's avatar
lambrechts committed
197
  _rho = NULL;
198
199
  _rhoGrad = NULL;
  _byPart = true;
lambrechts's avatar
lambrechts committed
200
201
202
203
204
}

void dgConservationLawSW3dPressure::setup() {
  if (!_rho)
    Msg::Fatal("dgConservationLawSW3dPressure: rho function must be set");
205
206
  if (_byPart){
    _volumeTerm0[""] = NULL;
207
208
    _volumeTerm1[""] = newFunctorMember(*this, &dgConservationLawSW3dPressure::advection);
    _interfaceTerm0[""] = newFunctorMember(*this, &dgConservationLawSW3dPressure::riemann);
209
210
  }
  else{
211
212
213
    if (!_rho)
      Msg::Fatal("dgConservationLawSW3dPressure: rho_grad function must be set");
    _volumeTerm0[""] = newFunctorMember(*this, &dgConservationLawSW3dPressure::source);
214
    _volumeTerm1[""] = NULL;
215
    _interfaceTerm0[""] = new functionConstant(std::vector<double>(2,0));
216
  }
lambrechts's avatar
lambrechts committed
217
218
219
220
221
222
223
}
dgConservationLawSW3dPressure::~dgConservationLawSW3dPressure() {
  if (_volumeTerm0[""]) delete _volumeTerm0[""];
  if (_volumeTerm1[""]) delete _volumeTerm1[""];
  if (_interfaceTerm0[""]) delete _interfaceTerm0[""];
}

224
/*class dgSW3dDepthIntegrate::source : public function {
lambrechts's avatar
lambrechts committed
225
226
227
  fullMatrix<double> _source;
  int _nbFields;
  public:
228
  source(const functor* source): function(source->getNbCol()) {
lambrechts's avatar
lambrechts committed
229
230
231
232
233
234
235
236
237
238
    setArgument (_source, source);
    _nbFields = source->getNbCol();
  };
  void call (dataCacheMap *m, fullMatrix<double> &val) {
    size_t nQP = val.size1();
    for(size_t i=0; i< nQP; i++) {
      for (int j=0; j < _nbFields; j++)
        val(i,j) = +_source(i,j);
    }
  }
239
};*/
lambrechts's avatar
lambrechts committed
240
241
242
243
244
245

class dgSW3dDepthIntegrate::riemann : public function {
  fullMatrix<double> normals,_solIn,_solOut;
  bool _propagateFromTop;
  int _nbFields;
  public:
246
  riemann(int nbF, const functor* source, bool propagateFromTop):function (2*nbF) {
lambrechts's avatar
lambrechts committed
247
248
249
250
    setArgument (normals, function::getNormals(), 0);
    setArgument (_solIn, function::getSolution(), 0);
    setArgument (_solOut, function::getSolution(), 1);
    _propagateFromTop=propagateFromTop;
251
    _nbFields = nbF;
lambrechts's avatar
lambrechts committed
252
253
254
255
  }
  void call (dataCacheMap *m, fullMatrix<double> &val) {
    // determine which element is below
    dataCacheMap* dgMap = m;
256
257
    int layerIdIn = dgMap->getGroupOfInterfaces()->elementGroup(0).getLayerId();
    int layerIdOut = dgMap->getGroupOfInterfaces()->elementGroup(1).getLayerId();
258
    bool outIsTop = layerIdIn > layerIdOut;
lambrechts's avatar
lambrechts committed
259
260
    bool isBoundary = dgMap->getGroupOfInterfaces()->nConnection() == 1;
    if (_propagateFromTop) // proper boundary behaviour
261
      outIsTop = layerIdIn >= layerIdOut;
lambrechts's avatar
lambrechts committed
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
    for(int i=0; i< val.size1(); i++) {
      double nz = normals(i,2);
      for (int j = 0; j < _nbFields; j++) {
        double solUpwindIn=0,solUpwindOut=0;
        if (!isBoundary) {
          if (outIsTop == _propagateFromTop) {
            solUpwindIn = -_solOut(i,j);
            solUpwindOut = 0;
          } else { //  out is top always
            solUpwindIn = 0;
            solUpwindOut = _solIn(i,j);
          }
        }else{
          if (layerIdIn == 0) {
            // bottom bnd condition
            if (_propagateFromTop) {
              solUpwindIn =  0;
              solUpwindOut = 0;
            } else {
              solUpwindIn = -_solOut(i,j);
              solUpwindOut = 0;
            }
          } else {
            // top bnd condition
            if (_propagateFromTop) {
              solUpwindIn =  -_solOut(i,j);
              solUpwindOut = 0;
            } else {
              solUpwindIn = 0;
              solUpwindOut = 0;
            }
          }
        }
        val(i,j) = solUpwindIn*nz;
        val(i,j+_nbFields) = solUpwindOut*nz;
      }
    }
  }
};

dgSW3dDepthIntegrate::dgSW3dDepthIntegrate( int nbFields, int propagateFromTop ):
  dgConservationLawFunction(nbFields), _nbFields(nbFields), _propagateFromTop(propagateFromTop)
{
  _f = NULL;
}

308
309
310
dgSW3dDepthIntegrate::dgSW3dDepthIntegrate(int nbF, const functor* func, int propagateFromTop):
  dgConservationLawFunction(nbF),
  _nbFields(nbF), _propagateFromTop(propagateFromTop),
lambrechts's avatar
lambrechts committed
311
312
313
314
315
316
317
318
319
  _f(func)
{
}

void dgSW3dDepthIntegrate::setup()
{
  if (!_f)
    Msg::Fatal("dgSW3dDepthIntegrate: Integrant function must be set");
    
320
  _volumeTerm0[""] = _f;
lambrechts's avatar
lambrechts committed
321
  _volumeTerm1[""] = NULL;
322
  _interfaceTerm0[""] = new riemann (_nbFields, _f,_propagateFromTop);
lambrechts's avatar
lambrechts committed
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
}

dgSW3dDepthIntegrate::~dgSW3dDepthIntegrate()
{
  if (_volumeTerm0[""]) delete _volumeTerm0[""];
  if (_volumeTerm1[""]) delete _volumeTerm1[""];
  if (_interfaceTerm0[""]) delete _interfaceTerm0[""];
}

dgFaceTermAssembler *dgSW3dDepthIntegrate::getFaceTerm(const dgGroupOfFaces &group) const
{
  if(group.orientationTag() == dgGroupOfFaces::ORIENTATION_HORIZONTAL)
    return dgConservationLawFunction::getFaceTerm(group);
  return NULL;
  
}


341

lambrechts's avatar
tmp    
lambrechts committed
342
343
344
dgConservationLawSW3dContinuitySolver::dgConservationLawSW3dContinuitySolver(dgDofContainer *sol, const dgExtrusion &extrusion, bool propagateFromTop):
  _extrusion(extrusion)
{
lambrechts's avatar
lambrechts committed
345
  _propagateFromTop = propagateFromTop;
346
  const dgGroupCollection &groups = *sol->getGroups();
lambrechts's avatar
lambrechts committed
347
348
  // build mass matrix
  fullMatrix<double> iMass;
349
350
351
352
353
354
355
356
  _iMass.resize(groups.getNbElementGroups());
  _neighbourRHS.resize(groups.getNbElementGroups());
  for(int iElemGroup=0;iElemGroup < groups.getNbElementGroups(); iElemGroup++){
    const dgGroupOfElements &group = *groups.getElementGroup(iElemGroup);
    const nodalBasis &fs = group.getFunctionSpace();
    int nbNodes = fs.points.size1();
    _iMass[iElemGroup].resize(nbNodes, nbNodes, group.getNbElements(), true);
    _neighbourRHS[iElemGroup].resize(group.getNbElements());
lambrechts's avatar
lambrechts committed
357
    const dgIntegrationMatrices &intMatrices = dgIntegrationMatrices::get(fs.type, fs.order * 2 + 1);
Sébastien Blaise's avatar
Sébastien Blaise committed
358
    const dgMeshJacobian &jacobians = groups.mesh().getJacobian(fs.order * 2 + 1);
lambrechts's avatar
lambrechts committed
359
360
361
    const fullMatrix<double> &detJ = jacobians.detJElement(iElemGroup);
    int nbQP = intMatrices.integrationPoints().size1();
    for (int iQP=0;iQP< nbQP ; iQP++ ){
362
      for (size_t iElem = 0; iElem < group.getNbElements(); iElem++){
lambrechts's avatar
lambrechts committed
363
364
365
        _iMass[iElemGroup].getBlockProxy(iElem, iMass);
        for (int k=0;k<nbNodes;k++){
          for (int l=0;l<nbNodes;l++) {
366
367
            for(int alpha = 0; alpha<group.getDimUVW(); alpha++)
              iMass(k,l) -= intMatrices.psiDPsiW()(k*nbNodes+l,alpha + group.getDimUVW()*iQP)*jacobians.invJElement(iElemGroup, iElem,iQP,alpha, 2)*detJ(iQP, iElem);
lambrechts's avatar
lambrechts committed
368
369
370
371
372
          }
        }
      }
    }
  }
373
374
375
376
  //add the vertical interface term that depends on the solution
  for(int iFaceGroup=0;iFaceGroup< groups.getNbFaceGroups();iFaceGroup++){
    const dgGroupOfFaces &faces = *groups.getFaceGroup(iFaceGroup);
    if (faces.orientationTag() != dgGroupOfFaces::ORIENTATION_HORIZONTAL)
lambrechts's avatar
lambrechts committed
377
      continue;
378
379
380
    const nodalBasis &fs = *faces.getNodalBasis();
    int nbFaceNodes = fs.points.size1();
    const dgIntegrationMatrices &intMatrices = dgIntegrationMatrices::get(fs.type, 2 * fs.order + 1);
Sébastien Blaise's avatar
Sébastien Blaise committed
381
    const dgMeshJacobian &jacobians = groups.mesh().getJacobian(2 * fs.order + 1);
382
    const fullMatrix<double> &detJ = jacobians.detJInterface(faces.interfaceVectorId());
lambrechts's avatar
lambrechts committed
383
384
    // propagateFromTop => take upper connection
    int nbQP = intMatrices.integrationPoints().size1();
385
386
387
388
389
390
391
    std::vector<int> iEGroup(faces.nConnection());
    for (int i = 0; i < faces.nConnection(); ++i)
      iEGroup[i] = groups.getElementGroupId(&faces.elementGroup(i));
    for (size_t iFace = 0 ; iFace < faces.size(); iFace++) {
      int iConn0 = 0, iConn1 = 0;
      if (faces.nConnection() == 2) {
        size_t col0, col1, lay0, lay1;
lambrechts's avatar
tmp    
lambrechts committed
392
393
        extrusion.mapElementID2ColumnPosition(iEGroup[0], faces.elementId(iFace, 0), col0, lay0);
        extrusion.mapElementID2ColumnPosition(iEGroup[1], faces.elementId(iFace, 1), col1, lay1);
394
395
396
        iConn0 = ((lay0 > lay1) == _propagateFromTop) ? 1 : 0;
        iConn1 = ((lay0 > lay1) == _propagateFromTop) ? 0 : 1;
      }
397
      const dgFullMatrix<double> &groupNormals = jacobians.normal(faces.interfaceVectorId(), iConn0);
398
399
400
401
402
403
      const std::vector<int> &closure0 = faces.closure(iFace, iConn0);
      const std::vector<int> &closure1 = faces.closure(iFace, iConn1);
      _iMass[iEGroup[iConn0]].getBlockProxy(faces.elementId(iFace, iConn0), iMass);
      fullMatrix<double> &rhsM = _neighbourRHS[iEGroup[iConn0]][faces.elementId(iFace, iConn0)];
      rhsM.resize(faces.elementGroup(iConn1).getNbNodes(), faces.elementGroup(iConn0).getNbNodes(), true);
      size_t faceLayer, dummy;
lambrechts's avatar
tmp    
lambrechts committed
404
      extrusion.mapDofToVEdgeDG(iEGroup[iConn0], faces.elementId(iFace, iConn0), closure0[0], dummy, faceLayer);
405
406
      if (faces.nConnection() == 1 && (_propagateFromTop == (faceLayer == 0)))
        continue;
lambrechts's avatar
lambrechts committed
407
      for (int iQP = 0; iQP < nbQP; iQP++) {
408
        double nz = groupNormals(iQP, 3 * faces.interfaceId(iFace) + 2);
lambrechts's avatar
lambrechts committed
409
410
        for (int k=0;k<nbFaceNodes;k++){
          for (int l=0;l<nbFaceNodes;l++) {
411
412
413
            double f =  nz * intMatrices.psiPsiW()(k*nbFaceNodes+l, iQP) * detJ(iQP, faces.interfaceId(iFace));
            iMass(closure0[l], closure0[k]) += f;
            rhsM(closure1[l], closure0[k]) += f;
lambrechts's avatar
lambrechts committed
414
415
416
417
418
419
420
421
          }
        }
      }
      iMass.invertInPlace();
    }
  }
}

422
void dgConservationLawSW3dContinuitySolver::iterate(dgConservationLawFunction *claw, dgDofContainer *sol, double time, bool scatterSolution){
423
  fullMatrix<double> residuEl, solEl, iMassEl;
lambrechts's avatar
lambrechts committed
424
  const dgGroupCollection &groups = *sol->getGroups();
425
426
427
428
429
430
  for (int i = 0; i < groups.getNbElementGroups(); ++i) {
    const dgGroupOfElements &group = *groups.getElementGroup(i);
    if (sol->nField(group) != claw->nField(group))
      Msg::Fatal("nbFields of the solution (%d) DofContainer does not match with the law (%d)",sol->nField(group),claw->nField(group));
  }
  dgDofContainer resd(*sol);
lambrechts's avatar
lambrechts committed
431
432
  // clear previous results (pollutes external values in Riemann)
  sol->setAll(0.0);
433
  claw->computeAllTerms(time, *sol, resd, false);
lambrechts's avatar
tmp    
lambrechts committed
434
435
436
  for (size_t iCol = 0; iCol < _extrusion.getNbColumns(); ++iCol) {
    for (size_t iL = 0; iL < _extrusion.getNbElemsInColumn(iCol); ++iL) {
      size_t iLay = _propagateFromTop ? iL : _extrusion.getNbElemsInColumn(iCol) - 1 - iL;
437
      size_t iGroup, iEl;
lambrechts's avatar
tmp    
lambrechts committed
438
      _extrusion.mapColumnPosition2ElementID(iCol, iLay, iGroup, iEl);
439
440
441
442
      resd.getGroupProxy(iGroup).getBlockProxy(iEl, residuEl);
      sol->getGroupProxy(iGroup).getBlockProxy(iEl, solEl);
      _iMass[iGroup].getBlockProxy(iEl, iMassEl);
      iMassEl.mult(residuEl, solEl);
lambrechts's avatar
tmp    
lambrechts committed
443
      if (iL != _extrusion.getNbElemsInColumn(iCol) - 1) {
444
445
        size_t iLay1 = _propagateFromTop ? iLay + 1 : iLay - 1;
        size_t iGroup1, iEl1;
lambrechts's avatar
tmp    
lambrechts committed
446
        _extrusion.mapColumnPosition2ElementID(iCol, iLay1, iGroup1, iEl1);
447
448
449
        resd.getGroupProxy(iGroup1).getBlockProxy(iEl1, residuEl);
        const fullMatrix<double > &M = _neighbourRHS[iGroup][iEl];
        residuEl.gemm(M, solEl);
lambrechts's avatar
lambrechts committed
450
451
452
453
454
455
456
      }
    }
  }
  if (scatterSolution)
    sol->scatter(true);
}

457
458
459
dgConservationLawSW3dContinuitySolver::~dgConservationLawSW3dContinuitySolver()
{
}
460
461
462
463
464

void dgConservationLawSW3dContinuity::finalize(const dgGroupOfElements &group, dgDofContainer &residual) const
{
  //dgConservationLawFunction::finalize(group, residual);
}
465
466
467
468
469

void dgConservationLawSW3dPressure::finalize(const dgGroupOfElements &group, dgDofContainer &residual) const
{
  //dgConservationLawFunction::finalize(group, residual);
}