dgDofContainer.cpp 95.5 KB
Newer Older
lambrechts's avatar
lambrechts committed
1
2
3
4
5
#include "dgDofContainer.h"
#include "dgGroupOfElements.h"
#include "dgIntegrationMatrices.h"
#include <limits.h>
#include <limits>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
6
#include "dgConfig.h"
lambrechts's avatar
lambrechts committed
7
8
9
10
11
12
13
14
15
16
17
18
#ifdef HAVE_MPI
#include "mpi.h"
#else
#include "string.h"
#endif
#include <sstream>
#include <fstream>
#include "dgDofManager.h"
#include "dgSpaceTransform.h"
#include "dgMesh.h"
#include "dgMeshJacobian.h"
#include "BasisFactory.h"
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
19
20
21
22
23
#include <algorithm>
#include "dgMessage.h"
#include "polynomialBasis.h"
#include "bezierBasis.h"
#include "OS.h"
Jonathan Lambrechts's avatar
HOmesh    
Jonathan Lambrechts committed
24
25
#include "SPoint3.h"
#include "SVector3.h"
26
#include "functor.h"
lambrechts's avatar
lambrechts committed
27

28
29
int dgDofContainer::nField(const dgGroupOfElements &group) const
{
30
31
32
  if (group.getGroupCollection() != _groups) {
    Msg::Fatal("groups does not match");
  }
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
  return _nField[group.elementVectorId()];
}

int dgDofContainer::nField(int elementVectorId) const
{
  return _nField[elementVectorId];
}

int dgDofContainer::nFieldMax() const
{
  return *std::max_element(_nField.begin(), _nField.end());
}

int dgDofContainer::nFieldConstant() const
{
  for (size_t i = 1; i < _nField.size(); ++i)
    if (_nField[i] != _nField[0])
      Msg::Fatal("this function requires a constant number of fields");
  return _nField[0];
}

lambrechts's avatar
lambrechts committed
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
class dgDofContainer::parallelStructure {
  public:
  class commChannel {
    public:
    int rank;
    std::vector<double> sendBuf, recvBuf;
    #ifdef HAVE_MPI
      MPI_Request sendRequest, recvRequest;
    #endif
    commChannel(int rank_, int nbSend, int nbRecv) {
      rank = rank_;
      sendBuf.resize(nbSend);
      recvBuf.resize(nbRecv);
    }
    void irecv(int tag)
    {
    #ifdef HAVE_MPI
      if (!recvBuf.empty())
        MPI_Irecv(&recvBuf[0], recvBuf.size(), MPI_DOUBLE, rank, tag, MPI_COMM_WORLD, &recvRequest);
    #endif
    }
    void isend(int tag)
    {
    #ifdef HAVE_MPI
      if (!sendBuf.empty())
        MPI_Isend(&sendBuf[0], sendBuf.size(), MPI_DOUBLE, rank, tag, MPI_COMM_WORLD, &sendRequest);
    #endif
    }
    void waitRecv()
    {
    #ifdef HAVE_MPI
      if (!recvBuf.empty())
        MPI_Wait(&recvRequest, MPI_STATUS_IGNORE);
    #endif
    }
    void waitSend()
    {
    #ifdef HAVE_MPI
      if (!sendBuf.empty())
        MPI_Wait(&sendRequest, MPI_STATUS_IGNORE);
    #endif
    }
  };
  int _maxTag;
  int _minTag;
  bool _scattering;
  bool _isDGGhost;
  std::map<int, commChannel> channels;
  dgDofContainer *_dof;

  parallelStructure(dgDofContainer *dof, bool isDGGhost, int maxTag = -1, int minTag = 0)
  {
    _scattering=false;
    _maxTag = maxTag;
    _minTag = minTag;
    _dof = dof;
    _isDGGhost = isDGGhost;
    const dgGroupCollection *groups = dof->getGroups();
    std::vector<size_t> nbSend(Msg::GetCommSize()), nbRecv(Msg::GetCommSize());
    for(int i=0;i<Msg::GetCommSize();i++) {
      nbSend[i] = 0;
      for(size_t j = 0; j < (size_t) groups->getNbImageElementsOnPartition(i, !_isDGGhost); j++){
        dgGroupOfElements *group= groups->getElementGroup(groups->getImageElementGroup(i,j, !_isDGGhost));
117
        int nbFields = dof->nField(*group);
lambrechts's avatar
lambrechts committed
118
119
120
121
122
123
        if (_maxTag == -1 || (_maxTag >= group->getMultirateTag() && group->getMultirateTag() >= _minTag))
          nbSend[i] += groups->getElementGroup(groups->getImageElementGroup(i,j, !_isDGGhost))->getNbNodes()*nbFields;
      }
      nbRecv[i] = 0;
      for(size_t j=0;j < (size_t) groups->getNbGhostElementsFromPartition(i, !_isDGGhost); j++){
        dgGroupOfElements *group= groups->getElementGroup(groups->getGhostElementGroup(i,j, !_isDGGhost));
124
        int nbFields = dof->nField(*group);
lambrechts's avatar
lambrechts committed
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
        if (_maxTag == -1 || (_maxTag >= group->getMultirateTag() && group->getMultirateTag() >= _minTag))
          nbRecv[i] += groups->getElementGroup(groups->getGhostElementGroup(i,j, !_isDGGhost))->getNbNodes() * nbFields;
      }
      if (nbSend[i] + nbRecv[i] > 0) {
        channels.insert(std::make_pair(i,commChannel(i, nbSend[i], nbRecv[i])));
      }
    }
  }
  void scatterBegin(int tag)
  {
    if (_scattering)
      Msg::Fatal("A scatter is already in progress on this dgDofContainer");
    _scattering=true;
    const dgGroupCollection *groups = _dof->getGroups();
    for (std::map<int,commChannel>::iterator it = channels.begin(); it != channels.end(); ++it)
      it->second.irecv(tag);
    for (std::map<int,commChannel>::iterator it = channels.begin(); it != channels.end(); ++it) {
      commChannel &channel = it->second;
      int index = 0;
      int rank = channel.rank;
      for(size_t j = 0; j < (size_t) groups->getNbImageElementsOnPartition(rank, !_isDGGhost); ++j){
146
147
148
149
        int groupId = groups->getImageElementGroup(rank, j, !_isDGGhost);
        const dgGroupOfElements &group = *groups->getElementGroup(groupId);
        int tag = group.getMultirateTag();
        int nbFields = _dof->nField(group);
lambrechts's avatar
lambrechts committed
150
151
        if (_maxTag != -1 && (_minTag > tag ||  tag > _maxTag))
          continue;
152
        fullMatrix<double> &sol = _dof->getGroupProxy(groupId);
lambrechts's avatar
lambrechts committed
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
        int eid = groups->getImageElementPositionInGroup(rank, j, !_isDGGhost);
        for(int l = 0; l < nbFields; ++l)
          for(int k = 0; k < sol.size1(); ++k)
            channel.sendBuf[index++] = sol(k, eid * nbFields + l);
      }
      channel.isend(tag);
    }
  }
  void scatterEnd()
  {
    if (!_scattering)
      Msg::Fatal("No scatter is in progress on this dgDofContainer");
    const dgGroupCollection *groups = _dof->getGroups();
    fullMatrix<double> solElProxy;
    for (std::map<int,commChannel>::iterator it = channels.begin(); it != channels.end(); ++it) {
      commChannel &channel = it->second;
      channel.waitRecv();
      int index = 0;
      int rank = channel.rank;
      for(size_t j = 0; j < (size_t) groups->getNbGhostElementsFromPartition(rank, !_isDGGhost); ++j){
lambrechts's avatar
lambrechts committed
173
        int groupId = groups->getGhostElementGroup(rank, j, !_isDGGhost);
174
175
176
        const dgGroupOfElements &group = *groups->getElementGroup(groupId);
        int tag = group.getMultirateTag();
        int nbFields = _dof->nField(group);
lambrechts's avatar
lambrechts committed
177
178
        if (_maxTag != -1 && (_minTag > tag ||  tag > _maxTag))
          continue;
179
        fullMatrix<double> &sol = _dof->getGroupProxy(groupId);
lambrechts's avatar
lambrechts committed
180
181
182
183
184
185
186
187
188
189
190
191
        int eid = groups->getGhostElementPositionInGroup(rank, j, !_isDGGhost);
        for(int l = 0; l < nbFields; ++l)
          for(int k = 0; k < sol.size1(); ++k)
            sol(k, eid * nbFields + l) = channel.recvBuf[index++];
      }
    }
    for (std::map<int,commChannel>::iterator it = channels.begin(); it != channels.end(); ++it)
      it->second.waitSend();
    _scattering=false;
  }
};

192
193


194
dgDofFunction::dgDofFunction(dgDofContainer *dof) :functor(){
195
196
  _dof = dof;
}
lambrechts's avatar
lambrechts committed
197

198
dgDofFunction::dgDofFunction():functor(), _dof(NULL){}
199

200
void dgDofFunction::setDof(dgDofContainer *dof) {
lambrechts's avatar
lambrechts committed
201
202
203
  _dof = dof;
}

204
205
206
207
208
209
210
211
212
213
214
215
static void transposeElement(const fullMatrix<double> &in, fullMatrix<double> &out, size_t nElement)
{
  out.resize(nElement * in.size1(), in.size2() / nElement, false);
  for (size_t i = 0; i < nElement; ++i) {
    for (int j = 0; j < in.size1(); ++j) {
      for (int k = 0; k < out.size2(); ++k) {
        out(i * in.size1() + j, k) = in(j, i * out.size2() + k);
      }
    }
  }
}

216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
//todo : we should store data in the right order on dgDofContainer...
static fullMatrix<double> transposeElementField(const dgFullMatrix<double> &in)
{
  fullMatrix<double> out(in.size1(), in.size2(), false);
  for (size_t iBlock = 0; iBlock < in.nBlock(); ++iBlock) {
    for (int j = 0; j < in.getBlockSize(); ++j) {
      for (int i = 0; i < in.size1(); ++i) {
        out(i, j * in.nBlock() + iBlock) = in(i, iBlock * in.getBlockSize() + j);
      }
    }
  }
  return out;
}


231
void dgDofFunction::operator()(functorCache &m, fullMatrix<double> &sol) const
lambrechts's avatar
lambrechts committed
232
{
233
234
235
236
  const int nFields = _dof->nField(*m.elementGroup());
  sol.resize(m.nEvaluationPoint(), nFields, false);
  switch(m.mode()) {
    case functorCache::INTEGRATION_GROUP_MODE:
lambrechts's avatar
lambrechts committed
237
    {
238
239
      if(m.interfaceGroup() == NULL) {
        _dof->_getAtIntegrationPoints(*m.elementGroup(), m.integrationMatrices(), sol);
lambrechts's avatar
lambrechts committed
240
      }
241
      else {
242
243
        if (m.side() >= m.interfaceGroup()->nConnection()) {
          m.invalidCall();
244
245
          return;
        }
246
247
        _dof->_getAtIntegrationPoints(*m.interfaceGroup(), m.side(),
                                     m.integrationMatrices(), sol);
lambrechts's avatar
lambrechts committed
248
249
250
      }
      return;
    }
251
    case functorCache::NODE_GROUP_MODE:
252
253
    {
      // nodes on element
254
255
      if(m.interfaceGroup() == NULL) {
        transposeElement(_dof->getGroupProxy(m.groupId()), sol, m.elementGroup()->getNbElements());
256
257
        return;
      }
258
259
      if (m.side() >= m.interfaceGroup()->nConnection()) {
        m.invalidCall();
260
261
        return;
      }
262
263

      // nodes on face
264
265
      const dgGroupOfFaces &faces = *m.interfaceGroup();
      const dgFullMatrix<double> &proxyGroup = _dof->getGroupProxy(m.groupId());
266
      for(size_t iF = 0; iF < faces.size(); ++iF) {
267
        const std::vector<int> &closure = m.interfaceGroup()->closure(m.interfaceId(iF), m.side());
268
        for (int iField = 0; iField < nFields; iField++) {
269
          const int ip = m.elementId(iF) * nFields + iField;
270
271
          for(size_t j = 0; j < closure.size(); j++) {
            sol(iF * closure.size() + j, iField) = proxyGroup(closure[j], ip);
272
273
274
275
276
          }
        }
      }
      return;
    }
277
    case functorCache::NODE_MODE:
lambrechts's avatar
lambrechts committed
278
279
    {
      // nodes on element
280
281
      if(m.interfaceGroup() == NULL) {
          _dof->getGroupProxy(m.groupId()).getBlockProxy(m.elementId(0), sol);
lambrechts's avatar
lambrechts committed
282
283
284
285
        return;
      }

      // nodes on face
286
287
      const dgFullMatrix<double> &proxyGroup = _dof->getGroupProxy(m.groupId());
      const std::vector<int> &closure = m.interfaceGroup()->closure(m.interfaceId(0), m.side());
288
      for (int iField = 0; iField < nFields; iField++) {
289
        const int ip = m.elementId(0) * nFields + iField;
290
291
        for(size_t j = 0; j < closure.size(); j++) {
          sol(j, iField) = proxyGroup(closure[j], ip);
lambrechts's avatar
lambrechts committed
292
293
294
295
296
        }
      }
      return;
    }

297
    case functorCache::POINT_MODE:
lambrechts's avatar
lambrechts committed
298
    {
299
      sol.setAll(0.);
lambrechts's avatar
lambrechts committed
300
301
      double f[256];
      fullMatrix<double> nodalV;
302
303
      const dgGroupOfFaces * faces = m.interfaceGroup();
      const fullMatrix<double> &uvw = m.parametricCoordinates();
lambrechts's avatar
lambrechts committed
304
305
      // volume case
      if(faces == NULL) {
306
        size_t groupId = m.groupId();
lambrechts's avatar
lambrechts committed
307
308
        const dgGroupOfElements *group = _dof->getGroups()->getElementGroup(groupId);
        group->getFunctionSpace().f(uvw(0, 0), uvw(0, 1), uvw(0, 2), f);
309
        _dof->getGroupProxy(groupId).getBlockProxy(m.elementId(0), nodalV);
lambrechts's avatar
lambrechts committed
310
        for (int j = 0; j < nodalV.size1(); j++) {
311
          for (int k = 0; k < nFields; k++) {
lambrechts's avatar
lambrechts committed
312
313
314
315
316
            sol(0, k) += nodalV(j, k) * f[j];
          }
        }
        return;
      }
317
      if (&m.groupCollection() != _dof->getGroups())
318
        Msg::Fatal("group collection does not match");
lambrechts's avatar
lambrechts committed
319
320
321
      // interface (boundary) case
      faces->getNodalBasis()->f(uvw(0, 0), uvw(0, 1), uvw(0, 2), f);
      int nbNodes = faces->getNbNodes();
322
323
      _dof->getGroupProxy(m.groupId()).getBlockProxy(m.elementId(0), nodalV);
      const std::vector<int> &closure = faces->closure(m.interfaceId(0), 0);
lambrechts's avatar
lambrechts committed
324
      for (int j = 0; j < nbNodes; j++) {
325
        for (int k = 0; k < nFields; k++) {
lambrechts's avatar
lambrechts committed
326
327
328
329
330
          sol(0, k) += nodalV(closure[j], k) * f[j];
        }
      }
      return;
    }
331
    case functorCache::MULTIPOINT_MODE:
332
333
334
335
    {
      sol.setAll(0.);
      double f[256];
      fullMatrix<double> nodalV;
336
337
      const dgGroupOfFaces * faces = m.interfaceGroup();
      const fullMatrix<double> &uvw = m.parametricCoordinates();
338
339
      // volume case
      if(faces == NULL) {
340
341
        size_t groupId = m.groupId();
        for (int i = 0; i < m.nEvaluationPoint(); ++i) {
342
343
          const dgGroupOfElements *group = _dof->getGroups()->getElementGroup(groupId);
          group->getFunctionSpace().f(uvw(i, 0), uvw(i, 1), uvw(i, 2), f);
344
          _dof->getGroupProxy(groupId).getBlockProxy(m.elementId(i), nodalV);
345
          for (int j = 0; j < nodalV.size1(); j++) {
346
            for (int k = 0; k < nFields; k++) {
347
348
349
350
351
352
              sol(i, k) += nodalV(j, k) * f[j];
            }
          }
        }
        return;
      }
353
      if (&m.groupCollection() != _dof->getGroups())
354
        Msg::Fatal("group collection does not match");
355
      // interface (boundary) case
356
      for (int i = 0; i < m.nEvaluationPoint(); ++i) {
357
358
        faces->getNodalBasis()->f(uvw(i, 0), uvw(i, 1), uvw(i, 2), f);
        int nbNodes = faces->getNbNodes();
359
360
        _dof->getGroupProxy(m.groupId()).getBlockProxy(m.elementId(i), nodalV);
        const std::vector<int> &closure = faces->closure(m.interfaceId(i), 0);
361
        for (int j = 0; j < nbNodes; j++) {
362
          for (int k = 0; k < nFields; k++) {
363
364
365
366
367
368
            sol(i, k) += nodalV(closure[j], k) * f[j];
          }
        }
      }
      return;
    }
lambrechts's avatar
lambrechts committed
369
370
371
372
373
    default:
      Msg::Fatal("a function requires a dof function but the algorithm does not provide it");
  }
}

374
dgDofFunctionGradient::dgDofFunctionGradient(dgDofContainer* dof):
375
  _dof(dof) {}
lambrechts's avatar
lambrechts committed
376

377
dgDofFunctionGradient::dgDofFunctionGradient(): _dof(NULL) {}
lambrechts's avatar
lambrechts committed
378
379
380
381
382

void dgDofFunctionGradient::setDof(dgDofContainer* dof) {
  _dof = dof;
}

383
void dgDofFunctionGradient::operator()(functorCache &m, fullMatrix< double >& sol)const
lambrechts's avatar
lambrechts committed
384
{
385
  sol.resize(m.nEvaluationPoint(), _dof->nField(*m.elementGroup())*3, false);
386
  if (m.mode() == functorCache::POINT_MODE) {
lambrechts's avatar
lambrechts committed
387
388
    fullMatrix<double> dVdXi, dVdX, S, dVdXENode;
    double f[256];
389
390
391
392
    const dgGroupOfElements *group = m.elementGroup();
    const dgIntegrationMatrices &intMatrices = m.integrationMatrices();
    const dgMeshJacobian &jac = m.jacobian();
    const fullMatrix<double> &uvw = m.parametricCoordinates();
lambrechts's avatar
lambrechts committed
393
    group->getFunctionSpace().f(uvw(0, 0), uvw(0, 1), uvw(0, 2), f);
394
395
    int nField = _dof->nField(*group);
    fullMatrix<double> nodalV(group->getNbNodes(), nField);;
396
    _dof->getGroupProxy(group).getBlockProxy(m.elementId(0), nodalV);
397
398
399
400
    dVdXi.resize(group->getDimUVW() * intMatrices.integrationPoints().size1(), nField);
    dVdX.resize(intMatrices.integrationPoints().size1(), 3 * nField);
    S.resize(group->getNbNodes(), 3 * nField);
    dVdXENode.resize(group->getNbNodes(), 3 * nField);
lambrechts's avatar
lambrechts committed
401
    intMatrices.dPsi().mult(nodalV, dVdXi); // dVdXi = V * dPsidXi
402
    jac.multDXiDXElement(m.groupId(), m.elementId(0), dVdXi, dVdX); // dVdX = dVdXi * dXidX
lambrechts's avatar
lambrechts committed
403
404
    intMatrices.psiW().mult(dVdX, S); // S = < dVdX * Psi >
    intMatrices.massInverse().mult(S, dVdXENode); // dVdXENode = M^-1 * S
405
    sol.resize(1, nField * 3, true);
lambrechts's avatar
lambrechts committed
406
    for (int j = 0; j < group->getNbNodes(); j++) {
407
      for (int k = 0; k < nField; k++) {
lambrechts's avatar
lambrechts committed
408
409
410
411
412
413
414
        sol(0, k * 3 + 0) += dVdXENode(j, k * 3 + 0) * f[j];
        sol(0, k * 3 + 1) += dVdXENode(j, k * 3 + 1) * f[j];
        sol(0, k * 3 + 2) += dVdXENode(j, k * 3 + 2) * f[j];
      }
    }
    return;
  }
415
  if (m.mode() == functorCache::MULTIPOINT_MODE) {
416
417
418
    sol.setAll(0.);
    fullMatrix<double> dVdXi, dVdX, S, dVdXENode;
    double f[256];
419
420
421
422
    const dgGroupOfElements *group = m.elementGroup();
    const dgIntegrationMatrices &intMatrices = m.integrationMatrices();
    const dgMeshJacobian &jac = m.jacobian();
    const fullMatrix<double> &uvw = m.parametricCoordinates();
423
    int nField = _dof->nField(*group);
424
    for (int i = 0; i < m.nEvaluationPoint(); ++i) {
425
      group->getFunctionSpace().f(uvw(i, 0), uvw(i, 1), uvw(i, 2), f);
426
      fullMatrix<double> nodalV(group->getNbNodes(), nField);;
427
      _dof->getGroupProxy(group).getBlockProxy(m.elementId(i), nodalV);
428
429
430
431
      dVdXi.resize(group->getDimUVW() * intMatrices.integrationPoints().size1(), nField);
      dVdX.resize(intMatrices.integrationPoints().size1(), 3 * nField);
      S.resize(group->getNbNodes(), 3 * nField);
      dVdXENode.resize(group->getNbNodes(), 3 * nField);
432
      intMatrices.dPsi().mult(nodalV, dVdXi); // dVdXi = V * dPsidXi
433
      jac.multDXiDXElement(m.groupId(), m.elementId(i), dVdXi, dVdX); // dVdX = dVdXi * dXidX
434
435
436
      intMatrices.psiW().mult(dVdX, S); // S = < dVdX * Psi >
      intMatrices.massInverse().mult(S, dVdXENode); // dVdXENode = M^-1 * S
      for (int j = 0; j < group->getNbNodes(); j++) {
437
        for (int k = 0; k < nField; k++) {
438
439
440
441
442
443
444
445
446
          sol(i, k * 3 + 0) += dVdXENode(j, k * 3 + 0) * f[j];
          sol(i, k * 3 + 1) += dVdXENode(j, k * 3 + 1) * f[j];
          sol(i, k * 3 + 2) += dVdXENode(j, k * 3 + 2) * f[j];
        }
      }
    }
    return;
    
  }
447
  if (m.mode() == functorCache::INTEGRATION_GROUP_MODE) {
448
    dgFullMatrix<double> groupval;
449
450
    if(m.interfaceGroup() == NULL) {
      _dof->_getGradientAtIntegrationPoints(*m.elementGroup(), m.integrationMatrices(), groupval);
lambrechts's avatar
lambrechts committed
451
    }
452
    else {
453
454
      if (m.side() >= m.interfaceGroup()->nConnection()) {
        m.invalidCall();
455
456
        return;
      }
457
      _dof->_getGradientAtIntegrationPoints(*m.interfaceGroup(), m.side(), m.integrationMatrices(), groupval);
458
459
460
461
    }
    transposeElement(groupval, sol, groupval.nBlock());
    return;
  }
462
  int nField = _dof->nField(*m.elementGroup());
463
  if (m.mode() == functorCache::NODE_GROUP_MODE) {
464
    dgFullMatrix<double> groupValue;
465
    if (m.interfaceGroup() == NULL) {
lambrechts's avatar
lambrechts committed
466
      dgFullMatrix<double> gradAtQP;
467
468
469
470
471
      _dof->_getGradientAtIntegrationPoints(*m.elementGroup(), m.integrationMatrices(), gradAtQP);
      fullMatrix<double> integratedGrad(m.elementGroup()->getNbNodes(), gradAtQP.size2());
      m.integrationMatrices().psiW().mult(gradAtQP, integratedGrad);
      groupValue.resize(integratedGrad.size1(), nField * 3, m.elementGroup()->getNbElements(), false);
      m.integrationMatrices().massInverse().mult(integratedGrad, groupValue);
lambrechts's avatar
lambrechts committed
472
    }
473
    else {
lambrechts's avatar
lambrechts committed
474
      dgFullMatrix<double> gradAtQP;
475
476
477
478
479
      _dof->_getGradientAtIntegrationPoints(*m.interfaceGroup(), m.side(), m.integrationMatrices(), gradAtQP);
      fullMatrix<double> integratedGrad(m.interfaceGroup()->getNodalBasis()->points.size1(), gradAtQP.size2());
      m.integrationMatrices().psiW().mult(gradAtQP, integratedGrad);
      groupValue.resize(integratedGrad.size1(), nField * 3, m.interfaceGroup()->size(), false);
      m.integrationMatrices().massInverse().mult(integratedGrad, groupValue);
lambrechts's avatar
lambrechts committed
480
    }
481
    transposeElement(groupValue, sol, groupValue.nBlock());
lambrechts's avatar
lambrechts committed
482
  }
483
  if(m.mode() == functorCache::NODE_MODE) {
484
485
    Msg::Fatal("node mode no longer supported");
    /*if (d->_groupValue.size1() == 0) {
486
      if (m->getGroupOfInterfaces() == NULL) {
487
        dgFullMatrix<double> gradAtQP;
488
489
490
491
492
        _dof->_getGradientAtIntegrationPoints(*m->getGroupOfElements(), m->getIntegrationMatrices(), gradAtQP);
        fullMatrix<double> integratedGrad(m->getGroupOfElements()->getNbNodes(), gradAtQP.size2());
        m->getIntegrationMatrices().psiW().mult(gradAtQP, integratedGrad);
        d->_groupValue.resize(integratedGrad.size1(), nField * 3, m->getGroupOfElements()->getNbElements(), false);
        m->getIntegrationMatrices().massInverse().mult(integratedGrad, d->_groupValue);
493
494
495
      }
      else {
        dgFullMatrix<double> gradAtQP;
496
497
498
499
500
        _dof->_getGradientAtIntegrationPoints(*m->getGroupOfInterfaces(), m->connectionId(), m->getIntegrationMatrices(), gradAtQP);
        fullMatrix<double> integratedGrad(m->getGroupOfInterfaces()->getNodalBasis()->points.size1(), gradAtQP.size2());
        m->getIntegrationMatrices().psiW().mult(gradAtQP, integratedGrad);
        d->_groupValue.resize(integratedGrad.size1(), nField * 3, m->getGroupOfInterfaces()->size(), false);
        m->getIntegrationMatrices().massInverse().mult(integratedGrad, d->_groupValue);
501
502
      }
    }
503
    d->_groupValue.getBlockProxy(m->getGroupOfInterfaces() == NULL ? m->elementId(0) : m->interfaceId(0), sol);*/
504
  }
lambrechts's avatar
lambrechts committed
505
506
507
508
509
510
}

void dgDofContainer::build(const dgGroupCollection *groups)
{
  int dataSize = 0;
  for (int i=0; i< groups->getNbElementGroups();i++){
511
512
    const dgGroupOfElements &group = *groups->getElementGroup(i);
    dataSize += group.getNbNodes() * group.getNbElements() * nField(group.elementVectorId());
lambrechts's avatar
lambrechts committed
513
514
515
  }

  // allocate the big vectors
516
  _data.resize(dataSize, false);
lambrechts's avatar
lambrechts committed
517
518
519
520
521
  // create proxys for each group
  int offset = 0;
  _dataProxys.resize(groups->getNbElementGroups()+groups->getNbGhostGroups());
  _groupFirstDofId.resize(groups->getNbElementGroups()+groups->getNbGhostGroups());
  for (int i=0;i<groups->getNbElementGroups();i++){
522
    const dgGroupOfElements *group = groups->getElementGroup(i);
lambrechts's avatar
lambrechts committed
523
524
    int nbNodes    = group->getNbNodes();
    int nbElements = group->getNbElements();
525
    _dataProxys[i].setAsProxy(&_data(offset), nbNodes, nField(*group), nbElements);
lambrechts's avatar
lambrechts committed
526
    _groupFirstDofId[i] = offset;
527
    offset += nbNodes * nField(*group) * nbElements;
lambrechts's avatar
lambrechts committed
528
529
530
531
532
533
  }

  //ghosts
  int totalNbElementsGhost = 0;
  int dataSizeGhost = 0;
  for (int i=0; i<groups->getNbGhostGroups(); i++){
534
535
536
    const dgGroupOfElements &group = *groups->getGhostGroup(i);
    totalNbElementsGhost += group.getNbElements();
    dataSizeGhost += group.getNbNodes() * nField(group) * group.getNbElements();
lambrechts's avatar
lambrechts committed
537
538
  }

539
  _ghostData.resize(dataSizeGhost, false);
lambrechts's avatar
lambrechts committed
540
541
542
543
544
545
  int ghostOffset=0;
  for (int i=0;i<groups->getNbGhostGroups();i++){
    dgGroupOfElements *group = groups->getGhostGroup(i);
    int nbNodes    = group->getNbNodes();
    int nbElements = group->getNbElements();
    int gid = i+groups->getNbElementGroups();
546
    _dataProxys[gid].setAsProxy(&_ghostData(ghostOffset), nbNodes, nField(*group), nbElements);
lambrechts's avatar
lambrechts committed
547
    _groupFirstDofId[gid] = offset;
548
549
    offset += nbNodes * nField(*group) * nbElements;
    ghostOffset += nbNodes * nField(*group) * nbElements;
lambrechts's avatar
lambrechts committed
550
551
552
553
554
555
556
557
558
  }
}

void dgDofContainer::_init() {
  _parallel = NULL;
  _parallelDG = NULL;
  _scatterOp = SCATTER_NONE;
  _function = new dgDofFunction(this);
  _functionGradient = new dgDofFunctionGradient(this);
559
  _fieldsName = new std::string[nFieldMax()];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
560
  _fs = &_groups->mesh().elementVector(0).functionSpace();
lambrechts's avatar
lambrechts committed
561
562
563
564
565
  build(_groups);
}

dgDofContainer::dgDofContainer (const dgDofContainer &orig):
  _groups (orig.getGroups()),
566
  _nField(orig._nField)
lambrechts's avatar
lambrechts committed
567
568
569
570
571
{
  _init();
}

dgDofContainer::dgDofContainer (const dgGroupCollection &groups,  int nbFields):
572
  _groups (&groups)
lambrechts's avatar
lambrechts committed
573
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
574
  _nField.resize(_groups->mesh().nElementVector() + _groups->mesh().nGhostElementVector(), nbFields);
575
576
577
  _init();
}

578
579
580
dgDofContainer::dgDofContainer(const dgGroupCollection &groups, const std::map<std::string, int> &nbFieldsForTag)
{
  _groups = &groups;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
581
  _nField.resize(groups.mesh().nElementVector() + _groups->mesh().nGhostElementVector(), 0);
582
  for (size_t i = 0; i < _nField.size(); ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
583
    const std::string &tag = groups.mesh().elementVector(i).physicalTag();
henrotte's avatar
henrotte committed
584
    std::map<std::string, int>::const_iterator it = nbFieldsForTag.find(tag);
585
586
587
588
589
590
591
592
    if (it == nbFieldsForTag.end())
      Msg::Warning("number of fields no specified for tag '%s', using 0 by default", tag.c_str());
    else
      _nField[i] = it->second;
  }
  _init();
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
593

594
595
596
597
#include "dgConservationLaw.h"
dgDofContainer::dgDofContainer (const dgGroupCollection &groups, const dgConservationLaw &law):
  _groups (&groups)
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
598
  _nField.resize(_groups->mesh().nElementVector() + _groups->mesh().nGhostElementVector());
599
  for (int i = 0; i < _groups->getNbElementGroups() + _groups->getNbGhostGroups(); ++i) {
600
601
602
    const dgGroupOfElements &group = * _groups->getElementGroup(i);
    _nField[group.elementVectorId()] = law.nField(group);
  }
lambrechts's avatar
lambrechts committed
603
604
  _init();
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
605

lambrechts's avatar
lambrechts committed
606
607
608
609
610
611
612

void dgDofContainer::transferGroups(const dgGroupCollection *newGroups, dgDofContainer *newDof){
  std::map<int, int> dataMap;
  int oldEls=0;
  for(int iGroup=0; iGroup<_groups->getNbElementGroups();iGroup++){
    dgGroupOfElements *g=_groups->getElementGroup(iGroup);
    for(size_t iElem=0; iElem<g->getNbElements();iElem++){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
613
      dataMap[g->getElement(iElem)->num()]=oldEls+iElem;
lambrechts's avatar
lambrechts committed
614
615
616
617
618
619
620
    }
    oldEls+=g->getNbElements();
  }

  int newEls=0;
  for(int iGroup=0; iGroup<newGroups->getNbElementGroups();iGroup++){
    dgGroupOfElements *g=newGroups->getElementGroup(iGroup);
621
    int nF = nField(*g);
lambrechts's avatar
lambrechts committed
622
623
624
    for(size_t iElem=0; iElem<g->getNbElements();iElem++){
      int nbNodes=g->getNbNodes();
      for(int iNode=0; iNode<nbNodes; iNode++){
625
        for(int iField=0; iField < nField(*g); iField++){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
626
          newDof->_data((newEls+iElem)*nbNodes*nF+iNode*nF+iField)=_data(dataMap[g->getElement(iElem)->num()]*nbNodes*nF+iNode*nF+iField);
lambrechts's avatar
lambrechts committed
627
628
629
630
631
632
633
634
        }
      }
    }
    newEls+=g->getNbElements();
  }
}

void dgDofContainer::changeGroups(const dgGroupCollection *newGroups, std::map<int, std::pair<int, int> > groupsMap){
635
636
637
  // the _mesh is not supposed to change but for now, it is changing, so we have to rebuild _nField
  std::map<std::string, int> str2nf;
  for (size_t i = 0; i < _nField.size(); ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
638
    str2nf[_groups->mesh().elementVector(i).physicalTag()] = _nField[i];
639
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
640
  _nField.resize(newGroups->mesh().nElementVector() + newGroups->mesh().nGhostElementVector(), 0);
641
  for (size_t i = 0; i < _nField.size(); ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
642
    _nField[i] = str2nf[newGroups->mesh().elementVector(i).physicalTag()];
643
  }
lambrechts's avatar
lambrechts committed
644
645
646
647
  fullVector<double> oldData=_data;
  build(newGroups);
  for(int iGroup=0; iGroup<newGroups->getNbElementGroups(); iGroup++){
    dgGroupOfElements *g=newGroups->getElementGroup(iGroup);
648
    int nF = nField(*g);
lambrechts's avatar
lambrechts committed
649
650
    int nbNodes=g->getNbNodes();
    for(size_t iElem=0; iElem<g->getNbElements();iElem++){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
651
652
      int oldJ=groupsMap[g->getElement(iElem)->num()].first;
      int newJ=groupsMap[g->getElement(iElem)->num()].second;
lambrechts's avatar
lambrechts committed
653
      for(int iNode=0; iNode<nbNodes; iNode++){
654
655
        for(int iField=0; iField<nF; iField++) {
          _data(newJ*nF+iNode*nF+iField) = oldData(oldJ*nF+iNode*nF+iField);
lambrechts's avatar
lambrechts committed
656
657
658
659
660
661
662
663
664
665
        }
      }
    }
  }
}

void dgDofContainer::setFieldName (int fieldIndex, std::string fieldName) {
  _fieldsName[fieldIndex] = fieldName;
}

vvallaeys's avatar
vvallaeys committed
666
667
668
669
std::string dgDofContainer::getFieldName (int fieldIndex) const {
  return _fieldsName[fieldIndex];
}

lambrechts's avatar
lambrechts committed
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
void dgDofContainer::destroyParallel(){
  if(_parallel){
    delete _parallel;
  }
  if(_parallelDG) {
    delete _parallelDG;
  }
  for (std::map<std::pair<int, int>, parallelStructure*>::iterator it = _parallelMRTag.begin(); it != _parallelMRTag.end(); ++it) {
    delete it->second;
  }
}

dgDofContainer::~dgDofContainer (){
  destroyParallel();
  delete _function;
  delete _functionGradient;
  delete[] _fieldsName;
}


void dgDofContainer::scatterBegin(bool dgGhost)
{
  static int ptag = 1000;
  ptag ++;
  if(_scatterOp != SCATTER_NONE) {
    Msg::Fatal("As scatter operation is already in progress");
  }
  if (dgGhost) {
    if (!_parallelDG)
      _parallelDG = new parallelStructure(this, true);
     _scatterOp = SCATTER_DG;
    _parallelDG->scatterBegin(ptag);
  }
  else {
    if (!_parallel)
      _parallel = new parallelStructure(this, false);
    _scatterOp = SCATTER_CG;
    _parallel->scatterBegin(ptag);
  }
}

void dgDofContainer::scatterBeginMRTag(int maxTag, int minTag)
{
  static int ptag = 1000;
  std::pair<int, int> tagPair = std::make_pair(maxTag, minTag);
  ptag ++;
  if (!_parallelMRTag[tagPair]){
    _parallelMRTag[tagPair] = new parallelStructure(this, true, maxTag, minTag);
  }
  _parallelMRTag[tagPair]->scatterBegin(ptag);
}

void dgDofContainer::scatterEndMRTag(int maxTag, int minTag)
{
  std::pair<int, int> tagPair = std::make_pair(maxTag, minTag);
  if (!_parallelMRTag[tagPair])
    Msg::Fatal ("no scatter begin for this tag");
  _parallelMRTag[tagPair]->scatterEnd();
}

void dgDofContainer::scatterEnd()
{
  switch (_scatterOp) {
    case SCATTER_CG :
      _parallel->scatterEnd();
      break;
    case SCATTER_DG :
      _parallelDG->scatterEnd();
      break;
    default:
      Msg::Warning("dgDofContainer::scatterEnd() called but no scatter operation in progress\n");
  }
  _scatterOp = SCATTER_NONE;
}

void dgDofContainer::scatter(bool dgGhost)
{
  scatterBegin(dgGhost);
  scatterEnd();
}

void dgDofContainer::setAll(double v) {
  for(int i=0;i<_data.size();i++)
    _data(i)=v;
  for(int i=0;i<_ghostData.size();i++)
    _ghostData(i)=v;
}
void dgDofContainer::setAll(std::vector<double> v) {
  fullMatrix<double> data, ghostData;
  for (int iGroup = 0; iGroup < _groups->getNbElementGroups(); iGroup++) {
    const dgGroupOfElements &group = *_groups->getElementGroup(iGroup);
761
762
763
    int nF = nField(group);
    if(nF != (int)v.size())
      Msg::Fatal("vector size must match nbfields in dgDofContainer setAll(vector)");
lambrechts's avatar
lambrechts committed
764
765
766
    for (size_t iElement = 0; iElement < group.getNbElements(); ++iElement) {
      getGroupProxy(iGroup).getBlockProxy(iElement, data);
      getGroupProxy(iGroup).getBlockProxy(iElement, ghostData);
767
      for(int k=0; k < nF; k++) {
lambrechts's avatar
lambrechts committed
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
        for(int i=0; i < group.getNbNodes(); i++) {
          data(i, k) = v[k];
          ghostData(i, k) = v[k];
        }
      }
    }
  }
}
void dgDofContainer::scale(double f)
{
  _data.scale(f);
  _ghostData.scale(f);
}

void dgDofContainer::multiply(const dgDofContainer &x, const dgDofContainer &y){
783
784
  fullMatrix<double> dataP, dataX, dataY;
  for (int iGroup = 0; iGroup < _groups->getNbElementGroups() + _groups->getNbGhostGroups(); iGroup++) {
lambrechts's avatar
lambrechts committed
785
    const dgGroupOfElements &group = *_groups->getElementGroup(iGroup);
786
787
788
    int nF = nField(group);
    if(nF != x.nField(group) || nF != y.nField(group))
      Msg::Fatal("vector size must match nbfields in dgDofContainer multiply(vector, vector)");
lambrechts's avatar
lambrechts committed
789
790
791
792
    for (size_t iElement = 0; iElement < group.getNbElements(); ++iElement) {
      getGroupProxy(iGroup).getBlockProxy(iElement, dataP);
      x.getGroupProxy(iGroup).getBlockProxy(iElement, dataX);
      y.getGroupProxy(iGroup).getBlockProxy(iElement, dataY);
793
      for(int k=0; k < nF; k++) {
lambrechts's avatar
lambrechts committed
794
795
796
797
798
799
800
801
        for(int i=0; i < group.getNbNodes(); i++) {
          dataP(i, k) = dataX(i, k) * dataY(i, k);
        }
      }
    }
  }
}

802
void dgDofContainer::divide(const dgDofContainer &x, const dgDofContainer &y, int iFieldY){
803
804
  fullMatrix<double> dataP, dataX, dataY;
  for (int iGroup = 0; iGroup < _groups->getNbElementGroups() + _groups->getNbGhostGroups(); iGroup++) {
805
    const dgGroupOfElements &group = *_groups->getElementGroup(iGroup);
806
807
808
    int nF = nField(group);
    if( iFieldY == -1 && (nF != x.nField(group) || nF != y.nField(group)))
      Msg::Fatal("vector size must match nbfields in dgDofContainer divide(vector, vector)");
809
810
811
812
    for (size_t iElement = 0; iElement < group.getNbElements(); ++iElement) {
      getGroupProxy(iGroup).getBlockProxy(iElement, dataP);
      x.getGroupProxy(iGroup).getBlockProxy(iElement, dataX);
      y.getGroupProxy(iGroup).getBlockProxy(iElement, dataY);
813
      for(int k=0; k < nF; k++) {
814
        for(int i=0; i < group.getNbNodes(); i++) {
815
816
817
818
819
820
          if (iFieldY == -1){
            dataP(i, k) = dataX(i, k) / dataY(i, k);
          }
          else{
            dataP(i, k) = dataX(i, k) / dataY(i, iFieldY);
          }
821
822
823
824
825
826
        }
      }
    }
  }
}

827
828
829
830
831
832
833
void dgDofContainer::fabs(){
  fullMatrix<double> data, ghostData;
  for (int iGroup = 0; iGroup < _groups->getNbElementGroups(); iGroup++) {
    const dgGroupOfElements &group = *_groups->getElementGroup(iGroup);
    for (size_t iElement = 0; iElement < group.getNbElements(); ++iElement) {
      getGroupProxy(iGroup).getBlockProxy(iElement, data);
      getGroupProxy(iGroup).getBlockProxy(iElement, ghostData);
834
      for(int k=0; k < nField(group); k++) {
835
836
837
838
839
840
841
842
843
        for(int i=0; i < group.getNbNodes(); i++) {
          data(i, k) = std::fabs(data(i, k));
          ghostData(i, k) = std::fabs(ghostData(i, k));
        }
      }
    }
  }
}

844
void dgDofContainer::dotProduct(const dgDofContainer &x, fullVector<double> &product){
lambrechts's avatar
lambrechts committed
845
  fullMatrix<double> dataP, dataX;
846
847
848
  if(nFieldMax() != x.nFieldMax())
    Msg::Fatal("vector size must match nbfields in dgdofcontainer dotproduct(vector)");
  fullVector<double> prod(nFieldMax());
lambrechts's avatar
lambrechts committed
849
  prod.setAll(0);
850
  product.resize(nFieldMax());
lambrechts's avatar
lambrechts committed
851
852
853
854
855
  for (int iGroup = 0; iGroup < _groups->getNbElementGroups(); iGroup++) {
    const dgGroupOfElements &group = *_groups->getElementGroup(iGroup);
    for (size_t iElement = 0; iElement < group.getNbElements(); ++iElement) {
      getGroupProxy(iGroup).getBlockProxy(iElement, dataP);
      x.getGroupProxy(iGroup).getBlockProxy(iElement, dataX);
856
      for(int k=0; k < nField(group); k++) {
lambrechts's avatar
lambrechts committed
857
        for(int i=0; i < group.getNbNodes(); i++) {
858
          prod(k) += dataP(i, k) * dataX(i, k);
lambrechts's avatar
lambrechts committed
859
860
861
862
863
        }
      }
    }
  }
  #if defined(HAVE_MPI)
864
  MPI_Allreduce(&prod(0), &product(0), nFieldMax(), MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
865
  #else
demaet's avatar
demaet committed
866
  product = prod;
lambrechts's avatar
lambrechts committed
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
  #endif
}


void dgDofContainer::axpy(const dgDofContainer &x, double a)
{
  _data.axpy(x._data,a);
  _ghostData.axpy(x._ghostData,a);
}

int dgDofContainer::norm0() const{
  int localNorm = _data.size();
#ifdef HAVE_MPI
  int globalNorm;
  MPI_Allreduce(&localNorm, &globalNorm, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
  return globalNorm;
#else
  return localNorm;
#endif
}

double dgDofContainer::norm() const{
#ifdef HAVE_MPI
  double localNorm = _data.norm();
  double localNorm2 = localNorm * localNorm;
  double globalNorm2;
  MPI_Allreduce(&localNorm2, &globalNorm2, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  double globalNorm = sqrt(globalNorm2);
  return globalNorm;
#else
  return _data.norm();
#endif
}

double dgDofContainer::norm(std::vector<dgGroupOfElements *>gV) {
  double norm=0;
  for(size_t i=0;i<gV.size();i++){
    dgGroupOfElements* g=gV[i];
    fullMatrix<double>&groupData=getGroupProxy(g);
    for(int r=0;r<groupData.size1();r++){
      for(int c=0;c<groupData.size1();c++){
        norm+=groupData(r,c)*groupData(r,c);
      }
    }
  }
  return norm;
}
void dgDofContainer::save(const std::string name) {
  FILE *f = fopen (name.c_str(),"wb");
  _data.binarySave(f);
  fclose(f);
}
void dgDofContainer::load(const std::string name) {
  FILE *f = fopen (name.c_str(),"rb");
  _data.binaryLoad(f);
  fclose(f);
}

static void multAddInverseMassMatrix ( const dgGroupOfElements & group, const dgFullMatrix<double> &residu, dgFullMatrix<double> &sol)
{
  fullMatrix<double> residuEl;
  fullMatrix<double> solEl;
  fullMatrix<double> iMassEl;
  for(size_t i = 0; i < group.getNbElements(); i++) {
    residu.getBlockProxy(i, residuEl);
    sol.getBlockProxy(i, solEl);
    group.getInverseMassMatrix().getBlockProxy(i, iMassEl);
    solEl.gemm(iMassEl, residuEl);
  }
}

938
void dgDofContainer::L2Projection(const functor *f, dgDofContainer *sol) {
939
  dgDofContainer rhs(*this);
940
  //rhs.scale(0.);
941
  functorCache cacheMap(functorCache::INTEGRATION_GROUP_MODE, _groups);
lambrechts's avatar
lambrechts committed
942
943
944
945
946
  if(sol)
    cacheMap.setSolutionFunction(sol->getFunction(), sol->getFunctionGradient());
  for (int iGroup = 0; iGroup < _groups->getNbElementGroups(); iGroup++) {
    const dgGroupOfElements &group = *_groups->getElementGroup(iGroup);
    cacheMap.setGroup(&group);
947
948
    const dgIntegrationMatrices &intMatrices = cacheMap.integrationMatrices();
    const fullMatrix<double> &val = cacheMap.get(*f);
949
    if (val.size2() == nField(group))
950
      rhs.getGroupProxy(iGroup).gemm(intMatrices.psiW(), cacheMap.jacobian().multJElement(group.elementVectorId(), val), 1., 0.);
951
952
953
954
955
956
957
    else {
      fullMatrix<double> tmp(val.size1(), nField(group), false);
      for (int i = 0; i < val.size1(); ++i) {
        for (int j = 0; j < nField(group); ++j) {
          tmp(i, j) = val(i, j);
        }
      }
958
      rhs.getGroupProxy(iGroup).gemm(intMatrices.psiW(), cacheMap.jacobian().multJElement(group.elementVectorId(), tmp), 1., 0.);
959
    }
lambrechts's avatar
lambrechts committed
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
    getGroupProxy(iGroup).setAll(0.);
    multAddInverseMassMatrix(group, rhs.getGroupProxy(iGroup), getGroupProxy(iGroup));
  }
}

L2ProjectionContinuous::L2ProjectionContinuous(dgDofManager* dofManager) :
      _dofManager(dofManager)
{
  _dofManager->fillSparsityPattern(false);

  // fill the left hand side with the mass matrix
  _dofManager->zeroMatrix();
  const dgGroupCollection &groups = _dofManager->getGroups();
  fullMatrix<double> massFull;
  for (int iGroup = 0; iGroup < groups.getNbElementGroups (); iGroup++) {
    dgGroupOfElements *group = groups.getElementGroup(iGroup);
    int nbNodes = group->getNbNodes();
977
978
    int nbFields = _dofManager->nField(*group);
    massFull.resize(nbNodes * nbFields, nbNodes * nbFields, true);
lambrechts's avatar
lambrechts committed
979
980
981
    fullMatrix<double> massE;
    for (size_t iElement = 0 ; iElement < group->getNbElements(); iElement++) {
      group->getMassMatrix().getBlockProxy(iElement, massE);
982
      for (int k = 0; k < nbFields; k++)
lambrechts's avatar
lambrechts committed
983
984
985
986
987
988
989
990
        for (int i = 0; i < nbNodes; i++)
          for (int j = 0; j < nbNodes; j++)
            massFull(k * nbNodes + i, k * nbNodes + j) = massE(i, j);
      _dofManager->assembleLHSMatrix(iGroup, iElement, iGroup, iElement, massFull);
    }
  }
}

991
void L2ProjectionContinuous::apply(dgDofContainer *output, const functor * f) {
lambrechts's avatar
lambrechts committed
992
993
994
995
996


  // fill the right hand side with the function f
  _dofManager->zeroRHS();
  const dgGroupCollection &groups = _dofManager->getGroups();
997
  functorCache cacheMap(functorCache::INTEGRATION_GROUP_MODE, &groups);
998
  dgFullMatrix<double> fJPsi;
lambrechts's avatar
lambrechts committed
999
1000
1001
1002
  fullMatrix<double> proxy;
  for (int iGroup = 0; iGroup < groups.getNbElementGroups (); iGroup++) {
    dgGroupOfElements *group = groups.getElementGroup(iGroup);
    cacheMap.setGroup(group);
1003
    const dgIntegrationMatrices &intMatrices = cacheMap.integrationMatrices();
1004
1005
1006
1007
    int nbFields = _dofManager->nField(*group);
    if(nbFields != output->nField(*group))
      Msg::Fatal("size of nbFields (%d) in L2ProjectionContinuous different than size of the dgDofContainer (%d) !", nbFields, output->nField(*group));
    fJPsi.resize(group->getNbNodes(), nbFields, group->getNbElements(), false);
1008
    fJPsi.gemm(intMatrices.psiW(), cacheMap.jacobian().multJElement(group->elementVectorId(), cacheMap.get(*f)), 1., 0.);
lambrechts's avatar
lambrechts committed
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
    for (size_t iElement = 0; iElement < group->getNbElements(); iElement++) {
      fJPsi.getBlockProxy(iElement, proxy);
      _dofManager->assembleRHSVector(iGroup, iElement, proxy);
    }
  }
  _dofManager->solve();
  // copy solution in the dof container
  fullMatrix<double> v, outputE;
  for (int iGroup = 0; iGroup < groups.getNbElementGroups (); iGroup++) {
    dgGroupOfElements *group = groups.getElementGroup(iGroup);
1019
    v.resize (group->getNbNodes() * output->nField(*group), 1);
lambrechts's avatar
lambrechts committed
1020
1021
1022
1023
1024
1025
1026
1027
1028
    for (size_t iElement = 0 ; iElement < group->getNbElements(); iElement++) {
      _dofManager->getSolution(iGroup, iElement, v);
      output->getGroupProxy(group).getBlockProxy(iElement, outputE);
      v.reshape(outputE.size1(), outputE.size2());
      outputE.setAll(v);
    }
  }
}

1029
void dgDofContainer::fromElsewhereL2Projection(const functor *f, const dgGroupCollection *gcEw, std::vector<std::string> tags, dgDofContainer *sol) {
lambrechts's avatar
lambrechts committed
1030
1031
1032
1033
1034
1035
1036
1037
  if(_groups->getDim() > gcEw->getDim())
    Msg::Error("could not project smth in n dimensions onto smth >n dimensions");
  bool sameDim = false;
  if(gcEw->getDim() == _groups->getDim())
    sameDim = true;
  else if(gcEw->getDim() - _groups->getDim() > 1)
    Msg::Error("could not project smth in n dimensions onto smth in smt different than n or n-1 dim");

1038
  dgDofContainer rhs(*this);
lambrechts's avatar
lambrechts committed
1039
  rhs.scale(0.);
1040
  functorCache cacheMapEw(functorCache::INTEGRATION_GROUP_MODE, gcEw);
lambrechts's avatar
lambrechts committed
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
  if(sol)
    cacheMapEw.setSolutionFunction(sol->getFunction(), sol->getFunctionGradient());
  fullMatrix<double> source;

  for (int iGroup = 0; iGroup < _groups->getNbElementGroups(); iGroup++) {
    const dgGroupOfElements &gLocal = *_groups->getElementGroup(iGroup);
    if(!tags.empty() && std::find(tags.begin(), tags.end(), gLocal.getPhysicalTag()) == tags.end())
      continue;

    int ig, ie;
    if(sameDim) {
      gcEw->find(gLocal.getElement(0), ig, ie);
      if(ig == -1)
        Msg::Error("elmt not found in dgVolumeOnSurfaceL2Projection !");
      const dgGroupOfElements * g = gcEw->getElementGroup(ig);
      cacheMapEw.setGroup(g);
    } else {
      gcEw->findFace(gLocal.getElement(0), ig, ie);
      if(ig == -1)
        Msg::Error("elmt not found in dgVolumeOnSurfaceL2Projection !");
      const dgGroupOfFaces * g = gcEw->getFaceGroup(ig);
      cacheMapEw.setInterfaceGroup(g);
    }

    const nodalBasis &fs = gLocal.getFunctionSpace();
1066
    const dgIntegrationMatrices &intMatrices = dgIntegrationMatrices::get(fs.type, cacheMapEw.integrationOrder());
lambrechts's avatar
lambrechts committed
1067
    int nbQP = intMatrices.integrationPoints().size1();
1068
    const dgMeshJacobian &jacobians = _groups->mesh().getJacobian(cacheMapEw.integrationOrder()); //TODO secure for != orders
lambrechts's avatar
lambrechts committed
1069

1070
    dgFullMatrix<double> Source(nbQP, nField(gLocal), gLocal.getNbElements(), false);
lambrechts's avatar
lambrechts committed
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083

    for (size_t iElmt = 0; iElmt < gLocal.getNbElements(); iElmt++) {
      if(sameDim) {
        gcEw->find(gLocal.getElement(iElmt), ig, ie);
        if(ie == -1)
          Msg::Error("elmt not found in dgVolumeOnSurfaceL2Projection !");
      } else {
        gcEw->findFace(gLocal.getElement(iElmt), ig, ie);
        if(ie == -1)
          Msg::Error("elmt not found in dgVolumeOnSurfaceL2Projection !");
      }

      Source.getBlockProxy(iElmt, source);
1084
      const fullMatrix<double> &sourceTermMatrix = cacheMapEw.get(*f);
lambrechts's avatar
lambrechts committed
1085
1086
      for (int iPt = 0; iPt < nbQP; ++iPt) {
        const double detJ = jacobians.detJElement(iGroup)(iPt, iElmt);
1087
        for (int k = 0; k < nField(gLocal); k++)
1088
          source(iPt, k) = sourceTermMatrix(ie *nbQP + iPt, k) * detJ;
lambrechts's avatar
lambrechts committed
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
      }
    }
    rhs.getGroupProxy(iGroup).gemm(intMatrices.psiW(), Source);
    getGroupProxy(iGroup).setAll(0.);
    multAddInverseMassMatrix(gLocal, rhs.getGroupProxy(iGroup), getGroupProxy(iGroup));
  }
}


void dgDofContainer::copy(const dgDofContainer& x)
{
1100
1101
  _data.setAll(x._data);
  _ghostData.setAll(x._ghostData);
lambrechts's avatar
lambrechts committed
1102
1103
1104
1105
1106
1107
1108
1109
1110
}

// multiplies the data array by the inverse mass matrix in place
void dgDofContainer::multiplyByInvMassMatrix()
{
  fullMatrix<double> dataEl, iMassEl;
  for(int k = 0; k < _groups->getNbElementGroups(); k++) {
    dgGroupOfElements *group = _groups->getElementGroup(k);
    int nbNodes = group->getNbNodes();
1111
    fullMatrix<double> tempElemData(nbNodes, nField(*group));
lambrechts's avatar
lambrechts committed
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
    for(size_t i = 0; i < group->getNbElements(); i++) {
      getGroupProxy(k).getBlockProxy(i, dataEl);
      group->getInverseMassMatrix().getBlockProxy(i, iMassEl);
      iMassEl.mult(dataEl,tempElemData);
      dataEl.setAll(tempElemData);
    }
  }
}

void dgDofContainer::multiplyByMassMatrix()
{
  fullMatrix<double> dataEl, massEl;
  for(int k = 0; k < _groups->getNbElementGroups(); k++) {
    dgGroupOfElements *group = _groups->getElementGroup(k);
    int nbNodes = group->getNbNodes();
1127
    fullMatrix<double> tempElemData(nbNodes, nField(*group));
lambrechts's avatar
lambrechts committed
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
    for(size_t i = 0; i < group->getNbElements(); i++) {
      getGroupProxy(k).getBlockProxy(i, dataEl);
      group->getMassMatrix().getBlockProxy(i, massEl);
      massEl.mult(dataEl, tempElemData);
      dataEl.setAll(tempElemData);
    }
  }
}

void dgDofContainer::multiplyByLocalMatrix(const dgGroupOfElements &group, const fullMatrix<double> mat)
{
  fullMatrix<double> dataEl;
  int nbNodes = group.getNbNodes();