VertexDoFMacroEdge.hpp 40.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/*
 * Copyright (c) 2017-2019 Daniel Drzisga, Dominik Thoennes, Marcus Mohr, Nils Kohl.
 *
 * This file is part of HyTeG
 * (see https://i10git.cs.fau.de/hyteg/hyteg).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
20
#pragma once
Dominik Thoennes's avatar
Dominik Thoennes committed
21

Dominik Thoennes's avatar
Dominik Thoennes committed
22
23
24
25
#include "hyteg/Levelinfo.hpp"
#include "hyteg/types/matrix.hpp"
#include "hyteg/p1functionspace/VertexDoFMemory.hpp"
#include "hyteg/p1functionspace/VertexDoFIndexing.hpp"
26

Dominik Thoennes's avatar
Dominik Thoennes committed
27
28
29
30
#include "hyteg/facedofspace/FaceDoFIndexing.hpp"
#include "hyteg/petsc/PETScWrapper.hpp"
#include "hyteg/indexing/Common.hpp"
#include "hyteg/primitives/Cell.hpp"
31
32
#include "hyteg/Algorithms.hpp"
#include "hyteg/indexing/DistanceCoordinateSystem.hpp"
33
#include "hyteg/sparseassembly/SparseMatrixProxy.hpp"
34
#include "hyteg/sparseassembly/VectorProxy.hpp"
Dominik Thoennes's avatar
Dominik Thoennes committed
35

36
#include "core/math/KahanSummation.h"
37
38
#include "core/DataTypes.h"

Marcus Mohr's avatar
Marcus Mohr committed
39
#ifdef DEBUG_ELEMENTWISE
Dominik Thoennes's avatar
Dominik Thoennes committed
40
#include "hyteg/format.hpp"
Marcus Mohr's avatar
Marcus Mohr committed
41
42
#endif

43
namespace hyteg {
44
45
namespace vertexdof {
namespace macroedge {
Dominik Thoennes's avatar
Dominik Thoennes committed
46

47
48
49
using walberla::real_c;
using indexing::Index;

50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
inline indexing::Index getIndexInNeighboringMacroCell( const indexing::Index&  vertexDoFIndexInMacroEdge,
                                                       const Edge&             edge,
                                                       const uint_t&           neighborCellID,
                                                       const PrimitiveStorage& storage,
                                                       const uint_t&           level )
{
   const Cell&  neighborCell = *( storage.getCell( edge.neighborCells().at( neighborCellID ) ) );
   const uint_t localEdgeID  = neighborCell.getLocalEdgeID( edge.getID() );

   const std::array< uint_t, 4 > localVertexIDsAtCell = algorithms::getMissingIntegersAscending< 2, 4 >(
       {neighborCell.getEdgeLocalVertexToCellLocalVertexMaps().at( localEdgeID ).at( 0 ),
        neighborCell.getEdgeLocalVertexToCellLocalVertexMaps().at( localEdgeID ).at( 1 )} );

   const auto indexInMacroCell = indexing::basisConversion(
       vertexDoFIndexInMacroEdge, localVertexIDsAtCell, {0, 1, 2, 3}, levelinfo::num_microvertices_per_edge( level ) );
   return indexInMacroCell;
}
67
68

inline Point3D coordinateFromIndex( const uint_t & level, const Edge & edge, const Index & index )
69
{
70
  const real_t  stepFrequency = 1.0 / levelinfo::num_microedges_per_edge( level );
71
72
73
74
  const Point3D step         = ( edge.getCoordinates()[1] - edge.getCoordinates()[0] ) * stepFrequency;
  return edge.getCoordinates()[0] + step * real_c( index.x() );
}

75
76
template<typename ValueType >
inline ValueType assembleLocal(const uint_t & level, uint_t pos, const Matrix3r& localMatrix,
77
78
                               double* src,
                               double* coeff,
Nils Kohl's avatar
Nils Kohl committed
79
                               const std::array< stencilDirection, 3 >& vertices,
Daniel Drzisga's avatar
Daniel Drzisga committed
80
                               const std::array<uint_t,3>& idx)
81
82
{

83
84
85
  ValueType meanCoeff = 1.0/3.0 * (coeff[vertexdof::macroedge::indexFromVertex( level, pos, vertices[0] )]
                                 + coeff[vertexdof::macroedge::indexFromVertex( level, pos, vertices[1] )]
                                 + coeff[vertexdof::macroedge::indexFromVertex( level, pos, vertices[2] )]);
86

Daniel Drzisga's avatar
Daniel Drzisga committed
87
  ValueType tmp;
88
89
90
  tmp  = localMatrix(idx[0],idx[0]) * src[vertexdof::macroedge::indexFromVertex( level, pos, vertices[0] )]
         + localMatrix(idx[0],idx[1]) * src[vertexdof::macroedge::indexFromVertex( level, pos, vertices[1] )]
         + localMatrix(idx[0],idx[2]) * src[vertexdof::macroedge::indexFromVertex( level, pos, vertices[2] )];
Daniel Drzisga's avatar
Daniel Drzisga committed
91
  return meanCoeff * tmp;
92
93
}

Daniel Drzisga's avatar
Daniel Drzisga committed
94
95
96
97
98
99
template<typename ValueType>
inline void assembleLocalStencil(uint_t level, uint_t pos, const Matrix3r& localMatrix,
                                 real_t* opr_data,
                                 real_t* coeff,
                                 const std::array< stencilDirection, 3 >& vertices,
                                 const std::array<uint_t,3>& idx)
100
101
{

Daniel Drzisga's avatar
Daniel Drzisga committed
102
103
104
  ValueType meanCoeff = 1.0/3.0 * (coeff[ vertexdof::macroedge::indexFromVertex( level, pos, vertices[ 0 ] ) ]
                                   + coeff[ vertexdof::macroedge::indexFromVertex( level, pos, vertices[ 1 ] ) ]
                                   + coeff[ vertexdof::macroedge::indexFromVertex( level, pos, vertices[ 2 ] ) ]);
105
106
107
108
109
110

  opr_data[vertexdof::stencilIndexFromVertex( vertices[0] )] += meanCoeff * localMatrix(idx[0],idx[0]);
  opr_data[vertexdof::stencilIndexFromVertex( vertices[1] )] += meanCoeff * localMatrix(idx[0],idx[1]);
  opr_data[vertexdof::stencilIndexFromVertex( vertices[2] )] += meanCoeff * localMatrix(idx[0],idx[2]);
}

111
112
113
114
115
116
117
118
119
120
121
122
123
template< typename ValueType >
inline void interpolate(const uint_t & level, Edge &edge,
                        const PrimitiveDataID< FunctionMemory< ValueType >, Edge> &edgeMemoryId,
                        const ValueType & scalar )
{
  size_t rowsize = levelinfo::num_microvertices_per_edge(level);
  auto edgeData =  edge.getData(edgeMemoryId)->getPointer( level );
  for (size_t i = 1; i < rowsize - 1; ++i)
  {
    edgeData[i] = scalar;
  }
}

124
125
template< typename ValueType >
inline void interpolate(const uint_t & level, Edge &edge,
126
                            const PrimitiveDataID< FunctionMemory< ValueType >, Edge> &edgeMemoryId,
127
                            const std::vector<PrimitiveDataID<FunctionMemory< ValueType >, Edge>> &srcIds,
128
                            const std::function<ValueType(const hyteg::Point3D &, const std::vector<ValueType>&)> &expr)
129
{
130
  ValueType * edgeData = edge.getData( edgeMemoryId )->getPointer( level );
131

132
133
134
  std::vector< ValueType * > srcPtr;
  for( const auto & src : srcIds )
  {
135
    srcPtr.push_back( edge.getData(src)->getPointer( level ) );
136
137
  }

138
  std::vector<ValueType> srcVector( srcIds.size() );
139

140
  Point3D xBlend;
141

142
  for ( const auto & it : vertexdof::macroedge::Iterator( level, 1 ) )
143
  {
144
145
    const Point3D coordinate = coordinateFromIndex( level, edge, it );
    const uint_t  idx        = vertexdof::macroedge::indexFromVertex( level, it.x(), stencilDirection::VERTEX_C );
146

147
148
149
    for ( uint_t k = 0; k < srcPtr.size(); ++k )
    {
      srcVector[ k ] = srcPtr[ k ][ idx ];
150
    }
151
152
    edge.getGeometryMap()->evalF(coordinate, xBlend);
    edgeData[ idx ] = expr( xBlend, srcVector );
Dominik Thoennes's avatar
Dominik Thoennes committed
153
154
155
  }
}

156

157
158
159
160
161
162
163
164
165
166
167
168

template< typename ValueType >
inline void swap( const uint_t & level, Edge & edge,
                  const PrimitiveDataID< FunctionMemory< ValueType >, Edge > & srcID,
                  const PrimitiveDataID< FunctionMemory< ValueType >, Edge > & dstID )
{
  auto srcData = edge.getData( srcID );
  auto dstData = edge.getData( dstID );
  srcData->swap( *dstData, level );
}


169
170
template< typename ValueType >
inline void assign( const uint_t & level, Edge &edge,
171
                   const std::vector<ValueType> &scalars,
172
173
                   const std::vector<PrimitiveDataID< FunctionMemory< ValueType >, Edge>> &srcIds,
                   const PrimitiveDataID< FunctionMemory< ValueType >, Edge> &dstId) {
174

175
  size_t rowsize = levelinfo::num_microvertices_per_edge(level);
Dominik Thoennes's avatar
Dominik Thoennes committed
176

Daniel Drzisga's avatar
Daniel Drzisga committed
177
  for (size_t i = 1; i < rowsize - 1; ++i) {
178
    ValueType tmp = scalars[0]*edge.getData(srcIds[0])->getPointer( level )[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )];
Dominik Thoennes's avatar
Dominik Thoennes committed
179

Daniel Drzisga's avatar
Daniel Drzisga committed
180
    for (size_t k = 1; k < srcIds.size(); ++k) {
181
      tmp += scalars[k]*edge.getData(srcIds[k])->getPointer( level )[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )];
Dominik Thoennes's avatar
Dominik Thoennes committed
182
183
    }

184
    edge.getData(dstId)->getPointer( level )[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] = tmp;
Dominik Thoennes's avatar
Dominik Thoennes committed
185
186
187
  }
}

188

189
190
191
192
193
194
195
196
197
198
199
200
201
202
template< typename ValueType >
inline void add(const uint_t & level,
                const Edge & edge,
                const ValueType & scalar,
                const PrimitiveDataID<FunctionMemory< ValueType >, Edge> & dstId )
{
  size_t rowsize = levelinfo::num_microvertices_per_edge(level);

  for (size_t i = 1; i < rowsize - 1; ++i)
  {
    edge.getData(dstId)->getPointer( level )[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] += scalar;
  }
}

203
204
template< typename ValueType >
inline void add( const uint_t & level, Edge &edge,
205
                const std::vector<ValueType> &scalars,
206
207
                const std::vector<PrimitiveDataID<FunctionMemory< ValueType >, Edge>> &srcIds,
                const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &dstId) {
208

209
  size_t rowsize = levelinfo::num_microvertices_per_edge(level);
Dominik Thoennes's avatar
Dominik Thoennes committed
210

Daniel Drzisga's avatar
Daniel Drzisga committed
211
  for (size_t i = 1; i < rowsize - 1; ++i) {
212
    auto tmp = ValueType( 0 );
Dominik Thoennes's avatar
Dominik Thoennes committed
213

Daniel Drzisga's avatar
Daniel Drzisga committed
214
    for (size_t k = 0; k < srcIds.size(); ++k) {
215
      tmp += scalars[k]*edge.getData(srcIds[k])->getPointer( level )[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )];
Dominik Thoennes's avatar
Dominik Thoennes committed
216
217
    }

218
    edge.getData(dstId)->getPointer( level )[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] += tmp;
Dominik Thoennes's avatar
Dominik Thoennes committed
219
220
221
  }
}

222

223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
template< typename ValueType >
inline void multElementwise( const uint_t & level, 
                             Edge &edge,
                             const std::vector<PrimitiveDataID<FunctionMemory< ValueType >, Edge>> &srcIds,
                             const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &dstId) {

  size_t rowsize = levelinfo::num_microvertices_per_edge(level);
  auto dst = edge.getData( dstId )->getPointer( level );

  for (size_t i = 1; i < rowsize - 1; ++i) 
  {
    const uint_t idx = vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C );
    ValueType tmp = edge.getData(srcIds[0])->getPointer( level )[idx];

    for (size_t k = 1; k < srcIds.size(); ++k) {
      tmp *= edge.getData(srcIds[k])->getPointer( level )[idx];
    }

    dst[idx] = tmp;
  }
}


246
template< typename ValueType >
247
inline ValueType dot( const uint_t & level, Edge &edge, const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &lhsMemoryId,
248
                  const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &rhsMemoryId) {
249

250
  walberla::math::KahanAccumulator< ValueType > scalarProduct;
251
  size_t rowsize = levelinfo::num_microvertices_per_edge(level);
Dominik Thoennes's avatar
Dominik Thoennes committed
252

Daniel Drzisga's avatar
Daniel Drzisga committed
253
  for (size_t i = 1; i < rowsize - 1; ++i) {
254
    scalarProduct += edge.getData(lhsMemoryId)->getPointer( level )[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )]
255
        * edge.getData(rhsMemoryId)->getPointer( level )[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )];
Dominik Thoennes's avatar
Dominik Thoennes committed
256
257
  }

258
  return scalarProduct.get();
Dominik Thoennes's avatar
Dominik Thoennes committed
259
260
}

261
template< typename ValueType >
262
inline ValueType sum( const uint_t & level, const Edge & edge, const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &dataID, const bool & absolute)
263
{
264
  auto sum = ValueType(0);
265
266
267
268
269
  size_t rowsize = levelinfo::num_microvertices_per_edge(level);

  auto data = edge.getData( dataID )->getPointer( level );

  for (size_t i = 1; i < rowsize - 1; ++i) {
270
271
272
273
274
275
276
277
    if ( absolute )
    {
      sum += std::abs( data[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] );
    }
    else
    {
      sum += data[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )];
    }
278
279
280
281
  }
  return sum;
}

282

283
284
template< typename ValueType >
inline void apply( const uint_t & level, Edge &edge, const PrimitiveDataID< StencilMemory< ValueType >, Edge> &operatorId,
285
286
                  const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &srcId,
                  const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &dstId, UpdateType update) {
287

Nils Kohl's avatar
Nils Kohl committed
288
  typedef stencilDirection sD;
289
  size_t rowsize = levelinfo::num_microvertices_per_edge(level);
Dominik Thoennes's avatar
Dominik Thoennes committed
290

291
292
293
  auto opr_data = edge.getData(operatorId)->getPointer( level );
  auto src = edge.getData(srcId)->getPointer( level );
  auto dst = edge.getData(dstId)->getPointer( level );
294

295
  ValueType tmp;
Dominik Thoennes's avatar
Dominik Thoennes committed
296

Nils Kohl's avatar
Nils Kohl committed
297
298
299
300
301
  for (size_t i = 1; i < rowsize - 1; ++i)
  {
    const auto stencilIdxW = vertexdof::macroedge::stencilIndexOnEdge( sD::VERTEX_W );
    const auto stencilIdxC = vertexdof::macroedge::stencilIndexOnEdge( sD::VERTEX_C );
    const auto stencilIdxE = vertexdof::macroedge::stencilIndexOnEdge( sD::VERTEX_E );
302

Nils Kohl's avatar
Nils Kohl committed
303
304
305
    const auto dofIdxW = vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_W );
    const auto dofIdxC = vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C );
    const auto dofIdxE = vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_E );
306

Nils Kohl's avatar
Nils Kohl committed
307
308
309
    tmp = opr_data[ stencilIdxW ] * src[ dofIdxW ]
      + opr_data[ stencilIdxC ] * src[ dofIdxC ]
      + opr_data[ stencilIdxE ] * src[ dofIdxE ];
310

Nils Kohl's avatar
Nils Kohl committed
311
    for ( uint_t neighborFace = 0; neighborFace < edge.getNumNeighborFaces(); neighborFace++ )
Nils Kohl's avatar
Nils Kohl committed
312
    {
Nils Kohl's avatar
Nils Kohl committed
313
314
315
316
317
318
319
      const auto stencilIdxWNeighborFace = vertexdof::macroedge::stencilIndexOnNeighborFace( sD::VERTEX_W, neighborFace );
      const auto stencilIdxENeighborFace = vertexdof::macroedge::stencilIndexOnNeighborFace( sD::VERTEX_E, neighborFace );
      const auto stencilWeightW = opr_data[ stencilIdxWNeighborFace ];
      const auto stencilWeightE = opr_data[ stencilIdxENeighborFace ];
      const auto dofIdxWNeighborFace = vertexdof::macroedge::indexFromVertexOnNeighborFace( level, i, neighborFace, sD::VERTEX_W );
      const auto dofIdxENeighborFace = vertexdof::macroedge::indexFromVertexOnNeighborFace( level, i, neighborFace, sD::VERTEX_E );
      tmp += stencilWeightW * src[dofIdxWNeighborFace] + stencilWeightE * src[dofIdxENeighborFace];
320
    }
Dominik Thoennes's avatar
Dominik Thoennes committed
321

Nils Kohl's avatar
Nils Kohl committed
322
    for ( uint_t neighborCell = 0; neighborCell < edge.getNumNeighborCells(); neighborCell++ )
Nils Kohl's avatar
Nils Kohl committed
323
    {
Nils Kohl's avatar
Nils Kohl committed
324
325
      const auto stencilIdx = vertexdof::macroedge::stencilIndexOnNeighborCell( neighborCell, edge.getNumNeighborFaces() );
      const auto dofIdx = vertexdof::macroedge::indexFromVertexOnNeighborCell( level, i, neighborCell, edge.getNumNeighborFaces() );
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
      tmp += opr_data[ stencilIdx ] * src[ dofIdx ];
    }

    if (update == Replace) {
      dst[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] = tmp;
    } else if (update == Add) {
      dst[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] += tmp;
    }
  }
}


template< typename ValueType >
inline void applyPointwise( const uint_t & level, const Edge &edge, const PrimitiveDataID< StencilMemory< ValueType >, Edge> &operatorId,
                  const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &srcId,
                  const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &dstId, UpdateType update) {

  typedef stencilDirection sD;
  size_t rowsize = levelinfo::num_microvertices_per_edge(level);

  auto opr_data = edge.getData(operatorId)->getPointer( level );
  auto src = edge.getData(srcId)->getPointer( level );
  auto dst = edge.getData(dstId)->getPointer( level );

  ValueType tmp;

  const auto stencilSize = 3 + 2 * edge.getNumNeighborFaces();

  for (size_t i = 1; i < rowsize - 1; ++i)
  {
    const auto offset = vertexdof::macroedge::innerIndex( level, i ) * stencilSize;
    const auto stencilIdxW = offset + vertexdof::macroedge::stencilIndexOnEdge( sD::VERTEX_W );
    const auto stencilIdxC = offset + vertexdof::macroedge::stencilIndexOnEdge( sD::VERTEX_C );
    const auto stencilIdxE = offset + vertexdof::macroedge::stencilIndexOnEdge( sD::VERTEX_E );

    const auto dofIdxW = vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_W );
    const auto dofIdxC = vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C );
    const auto dofIdxE = vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_E );

    tmp = opr_data[ stencilIdxW ] * src[ dofIdxW ]
      + opr_data[ stencilIdxC ] * src[ dofIdxC ]
      + opr_data[ stencilIdxE ] * src[ dofIdxE ];

    for ( uint_t neighborFace = 0; neighborFace < edge.getNumNeighborFaces(); neighborFace++ )
    {
      const auto stencilIdxWNeighborFace = offset + vertexdof::macroedge::stencilIndexOnNeighborFace( sD::VERTEX_W, neighborFace );
      const auto stencilIdxENeighborFace = offset + vertexdof::macroedge::stencilIndexOnNeighborFace( sD::VERTEX_E, neighborFace );
      const auto stencilWeightW = opr_data[ stencilIdxWNeighborFace ];
      const auto stencilWeightE = opr_data[ stencilIdxENeighborFace ];
      const auto dofIdxWNeighborFace = vertexdof::macroedge::indexFromVertexOnNeighborFace( level, i, neighborFace, sD::VERTEX_W );
      const auto dofIdxENeighborFace = vertexdof::macroedge::indexFromVertexOnNeighborFace( level, i, neighborFace, sD::VERTEX_E );
      tmp += stencilWeightW * src[dofIdxWNeighborFace] + stencilWeightE * src[dofIdxENeighborFace];
    }

    for ( uint_t neighborCell = 0; neighborCell < edge.getNumNeighborCells(); neighborCell++ )
    {
      const auto stencilIdx = vertexdof::macroedge::stencilIndexOnNeighborCell( neighborCell, edge.getNumNeighborFaces() );
      const auto dofIdx = vertexdof::macroedge::indexFromVertexOnNeighborCell( level, i, neighborCell, edge.getNumNeighborFaces() );
Nils Kohl's avatar
Nils Kohl committed
384
      tmp += opr_data[ stencilIdx ] * src[ dofIdx ];
385
386
387
    }

    if (update == Replace) {
388
      dst[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] = tmp;
389
    } else if (update == Add) {
390
      dst[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] += tmp;
Dominik Thoennes's avatar
Dominik Thoennes committed
391
392
393
394
    }
  }
}

395

396
397
template< typename ValueType >
inline void smooth_gs( const uint_t & level, Edge &edge, const PrimitiveDataID< StencilMemory< ValueType >, Edge> &operatorId,
398
399
                      const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &dstId,
                      const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &rhsId) {
400

401
  typedef stencilDirection sD;
402
  size_t rowsize = levelinfo::num_microvertices_per_edge(level);
Dominik Thoennes's avatar
Dominik Thoennes committed
403

404
405
  auto opr_data = edge.getData(operatorId)->getPointer( level );
  auto rhs = edge.getData(rhsId)->getPointer( level );
406
  auto dst = edge.getData(dstId)->getPointer( level );
Dominik Thoennes's avatar
Dominik Thoennes committed
407

408
409
410
  const auto stencilIdxW = vertexdof::macroedge::stencilIndexOnEdge( sD::VERTEX_W );
  const auto stencilIdxC = vertexdof::macroedge::stencilIndexOnEdge( sD::VERTEX_C );
  const auto stencilIdxE = vertexdof::macroedge::stencilIndexOnEdge( sD::VERTEX_E );
411

412
  const auto invCenterWeight = 1.0 / opr_data[ stencilIdxC ];
413

414
415
416
417
418
419
420
  ValueType tmp;

  for (size_t i = 1; i < rowsize - 1; ++i)
  {
    const auto dofIdxW = vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_W );
    const auto dofIdxC = vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C );
    const auto dofIdxE = vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_E );
421

422
423
424
425
426
    tmp = rhs[ dofIdxC ];

    tmp -= opr_data[ stencilIdxW ] * dst[ dofIdxW ] + opr_data[ stencilIdxE ] * dst[ dofIdxE ];

    for ( uint_t neighborFace = 0; neighborFace < edge.getNumNeighborFaces(); neighborFace++ )
Nils Kohl's avatar
Nils Kohl committed
427
    {
Nils Kohl's avatar
Nils Kohl committed
428
429
430
431
432
433
434
      const auto stencilIdxWNeighborFace = vertexdof::macroedge::stencilIndexOnNeighborFace( sD::VERTEX_W, neighborFace );
      const auto stencilIdxENeighborFace = vertexdof::macroedge::stencilIndexOnNeighborFace( sD::VERTEX_E, neighborFace );
      const auto stencilWeightW = opr_data[ stencilIdxWNeighborFace ];
      const auto stencilWeightE = opr_data[ stencilIdxENeighborFace ];
      const auto dofIdxWNeighborFace = vertexdof::macroedge::indexFromVertexOnNeighborFace( level, i, neighborFace, sD::VERTEX_W );
      const auto dofIdxENeighborFace = vertexdof::macroedge::indexFromVertexOnNeighborFace( level, i, neighborFace, sD::VERTEX_E );
      tmp -= stencilWeightW * dst[dofIdxWNeighborFace] + stencilWeightE * dst[dofIdxENeighborFace];
435
    }
Dominik Thoennes's avatar
Dominik Thoennes committed
436

437
    for ( uint_t neighborCell = 0; neighborCell < edge.getNumNeighborCells(); neighborCell++ )
Nils Kohl's avatar
Nils Kohl committed
438
    {
439
440
441
      const auto stencilIdx = vertexdof::macroedge::stencilIndexOnNeighborCell( neighborCell, edge.getNumNeighborFaces() );
      const auto dofIdx = vertexdof::macroedge::indexFromVertexOnNeighborCell( level, i, neighborCell, edge.getNumNeighborFaces() );
      tmp -= opr_data[ stencilIdx ] * dst[ dofIdx ];
Dominik Thoennes's avatar
Dominik Thoennes committed
442
443
    }

444
    dst[ dofIdxC ] = invCenterWeight * tmp;
Dominik Thoennes's avatar
Dominik Thoennes committed
445
446
447
  }
}

448
449
template< typename ValueType >
inline void smooth_sor(const uint_t & level, Edge &edge, const PrimitiveDataID< StencilMemory< ValueType >, Edge> &operatorId,
450
451
                          const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &dstId,
                          const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &rhsId,
452
453
454
                          ValueType relax,
                          const bool & backwards = false )
{
455

Nils Kohl's avatar
Nils Kohl committed
456
  typedef stencilDirection sD;
457
  size_t rowsize = levelinfo::num_microvertices_per_edge(level);
458

459
460
  auto opr_data = edge.getData(operatorId)->getPointer( level );
  auto rhs = edge.getData(rhsId)->getPointer( level );
Nils Kohl's avatar
Nils Kohl committed
461
462
463
464
465
466
467
  auto dst = edge.getData(dstId)->getPointer( level );

  const auto stencilIdxW = vertexdof::macroedge::stencilIndexOnEdge( sD::VERTEX_W );
  const auto stencilIdxC = vertexdof::macroedge::stencilIndexOnEdge( sD::VERTEX_C );
  const auto stencilIdxE = vertexdof::macroedge::stencilIndexOnEdge( sD::VERTEX_E );

  const auto invCenterWeight = 1.0 / opr_data[ stencilIdxC ];
468
469
470

  ValueType tmp;

471
472
473
474
475
  const int start = backwards ? (int)rowsize - 2 : 1;
  const int stop = backwards ? 0 : (int)rowsize - 1;
  const int incr = backwards ? -1 : 1;

  for (int ii = start; ii != stop; ii += incr)
Nils Kohl's avatar
Nils Kohl committed
476
  {
477
    const uint_t i = uint_c(ii);
Nils Kohl's avatar
Nils Kohl committed
478
479
480
    const auto dofIdxW = vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_W );
    const auto dofIdxC = vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C );
    const auto dofIdxE = vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_E );
481

Nils Kohl's avatar
Nils Kohl committed
482
    tmp = rhs[ dofIdxC ];
483

Nils Kohl's avatar
Nils Kohl committed
484
    tmp -= opr_data[ stencilIdxW ] * dst[ dofIdxW ] + opr_data[ stencilIdxE ] * dst[ dofIdxE ];
485

Nils Kohl's avatar
Nils Kohl committed
486
    for ( uint_t neighborFace = 0; neighborFace < edge.getNumNeighborFaces(); neighborFace++ )
Nils Kohl's avatar
Nils Kohl committed
487
    {
Nils Kohl's avatar
Nils Kohl committed
488
489
490
491
492
493
494
      const auto stencilIdxWNeighborFace = vertexdof::macroedge::stencilIndexOnNeighborFace( sD::VERTEX_W, neighborFace );
      const auto stencilIdxENeighborFace = vertexdof::macroedge::stencilIndexOnNeighborFace( sD::VERTEX_E, neighborFace );
      const auto stencilWeightW = opr_data[ stencilIdxWNeighborFace ];
      const auto stencilWeightE = opr_data[ stencilIdxENeighborFace ];
      const auto dofIdxWNeighborFace = vertexdof::macroedge::indexFromVertexOnNeighborFace( level, i, neighborFace, sD::VERTEX_W );
      const auto dofIdxENeighborFace = vertexdof::macroedge::indexFromVertexOnNeighborFace( level, i, neighborFace, sD::VERTEX_E );
      tmp -= stencilWeightW * dst[dofIdxWNeighborFace] + stencilWeightE * dst[dofIdxENeighborFace];
495
496
    }

Nils Kohl's avatar
Nils Kohl committed
497
    for ( uint_t neighborCell = 0; neighborCell < edge.getNumNeighborCells(); neighborCell++ )
Nils Kohl's avatar
Nils Kohl committed
498
    {
Nils Kohl's avatar
Nils Kohl committed
499
500
501
      const auto stencilIdx = vertexdof::macroedge::stencilIndexOnNeighborCell( neighborCell, edge.getNumNeighborFaces() );
      const auto dofIdx = vertexdof::macroedge::indexFromVertexOnNeighborCell( level, i, neighborCell, edge.getNumNeighborFaces() );
      tmp -= opr_data[ stencilIdx ] * dst[ dofIdx ];
502
    }
Dominik Thoennes's avatar
Dominik Thoennes committed
503

Nils Kohl's avatar
Nils Kohl committed
504
    dst[ dofIdxC ] = (1.0 - relax) * dst[ dofIdxC ] + relax * invCenterWeight * tmp;
505
506
507
508
  }
}


509
510
template< typename ValueType >
inline void smooth_jac(const uint_t & level, Edge &edge, const PrimitiveDataID< StencilMemory< ValueType >, Edge> &operatorId,
511
512
513
                          const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &dstId,
                          const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &rhsId,
                          const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &tmpId) {
514

515
  size_t rowsize = levelinfo::num_microvertices_per_edge(level);
516

517
518
519
520
  auto opr_data = edge.getData(operatorId)->getPointer( level );
  auto dst = edge.getData(dstId)->getPointer( level );
  auto rhs = edge.getData(rhsId)->getPointer( level );
  auto tmp = edge.getData(tmpId)->getPointer( level );
521
522
523

  for (size_t i = 1; i < rowsize - 1; ++i) {

524
    dst[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] = rhs[vertexdof::macroedge::indexFromVertex( level, i,
525
                                                                                                                                                                stencilDirection::VERTEX_C )];
526

Nils Kohl's avatar
Nils Kohl committed
527
528
    for (auto& neighbor : vertexdof::macroedge::neighborsOnEdgeFromVertexDoF )
    {
529
      dst[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] -= opr_data[ vertexdof::stencilIndexFromVertex( neighbor ) ] * tmp[vertexdof::macroedge::indexFromVertex( level, i, neighbor )];
530
531
    }

Nils Kohl's avatar
Nils Kohl committed
532
533
    for (auto& neighbor : vertexdof::macroedge::neighborsOnSouthFaceFromVertexDoF )
    {
534
      dst[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] -= opr_data[ vertexdof::stencilIndexFromVertex( neighbor ) ] * tmp[vertexdof::macroedge::indexFromVertex( level, i, neighbor )];
535
536
    }

Nils Kohl's avatar
Nils Kohl committed
537
538
539
540
    if (edge.getNumNeighborFaces() == 2)
    {
      for (auto& neighbor : vertexdof::macroedge::neighborsOnNorthFaceFromVertexDoF )
      {
541
        dst[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] -= opr_data[ vertexdof::stencilIndexFromVertex( neighbor ) ] * tmp[vertexdof::macroedge::indexFromVertex( level, i, neighbor )];
542
543
544
      }
    }

545
    dst[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] /= opr_data[vertexdof::stencilIndexFromVertex( stencilDirection::VERTEX_C ) ];
546
547
548
549
  }
}


550
template< typename ValueType >
551
inline void enumerate( const uint_t & level, Edge &edge, const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &dstId, ValueType& num) {
552

553
  size_t rowsize = levelinfo::num_microvertices_per_edge(level);
Daniel Drzisga's avatar
Daniel Drzisga committed
554
555

  for (size_t i = 1; i < rowsize - 1; ++i) {
556
557
    edge.getData(dstId)->getPointer( level )[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] = num;
    num++;
Daniel Drzisga's avatar
Daniel Drzisga committed
558
559
560
  }
}

561

562
563
template< typename ValueType >
inline void integrateDG(const uint_t & level, Edge &edge,
564
565
                            const std::shared_ptr< PrimitiveStorage >& storage,
                            const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &rhsId,
566
                            const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &rhsP1Id,
567
568
569
570
                            const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &dstId) {

  typedef stencilDirection sD;

571
  size_t rowsize = levelinfo::num_microvertices_per_edge(level);
572

573
574
575
  auto rhs = edge.getData(rhsId)->getPointer( level );
  auto rhsP1 = edge.getData(rhsP1Id)->getPointer( level );
  auto dst = edge.getData(dstId)->getPointer( level );
576

577
  real_t tmp;
578
579
580

  Face* face = storage->getFace(edge.neighborFaces()[0]);

581
582
  real_t weightedFaceArea0 = std::pow(4.0, -walberla::real_c(level)) * face->area / 3.0;
  real_t weightedFaceArea1 = real_c( 0 );
583
584
585

  if (edge.getNumNeighborFaces() == 2) {
    face = storage->getFace(edge.neighborFaces()[1]);
586
     weightedFaceArea1 = std::pow(4.0, -walberla::real_c(level)) * face->area / 3.0;
587
588
589
590
  }

  for (size_t i = 1; i < rowsize - 1; ++i) {

591
592
593
594
595
596
597
598
599
600
601
602
    tmp = weightedFaceArea0 * rhs[facedof::macroedge::indexFaceFromVertex( level, i, sD::CELL_GRAY_SW )] * ( 0.5 * 0.5 * ( rhsP1[vertexdof::macroedge::indexFromVertex( level, i,
                                                                                                                                                                        sD::VERTEX_C )]
                                                                                                                           + rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_W )] )
                                                                                                             + 0.5 * 0.5 * ( rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C )] +
                                                                                                                             rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_S )] ));
    tmp += weightedFaceArea0 * rhs[facedof::macroedge::indexFaceFromVertex( level, i, sD::CELL_BLUE_SE )] * ( 0.5 * 0.5 * ( rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C )] +
                                                                                                                            rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_S )] ) +
                                                                                                              0.5 * 0.5 * ( rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C )] +
                                                                                                                            rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_SE )] ));
    tmp += weightedFaceArea0 * rhs[facedof::macroedge::indexFaceFromVertex( level, i, sD::CELL_GRAY_SE )] *
      ( 0.5 * 0.5 * ( rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C )] + rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_SE )] ) +
        0.5 * 0.5 * ( rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C )] + rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_E )] ));
603
604
605

    if (edge.getNumNeighborFaces() == 2) {

606
607
608
609
610
611
612
613
614
615
616
617
      tmp += weightedFaceArea1 * rhs[facedof::macroedge::indexFaceFromVertex( level, i, sD::CELL_GRAY_NW )] * ( 0.5 * 0.5 * ( rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C )]
                                                                                                                              + rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_W )])
                                                                                                                + 0.5 * 0.5 * ( rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C )]
                                                                                                                                + rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_NW )]));
      tmp += weightedFaceArea1 * rhs[facedof::macroedge::indexFaceFromVertex( level, i, sD::CELL_BLUE_NW )] * ( 0.5 * 0.5 * ( rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C )]
                                                                                                                              + rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_NW )])
                                                                                                                + 0.5 * 0.5 * ( rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C )]
                                                                                                                                + rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_N )]));
      tmp += weightedFaceArea1 * rhs[facedof::macroedge::indexFaceFromVertex( level, i, sD::CELL_GRAY_NE )] * ( 0.5 * 0.5 * ( rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C )]
                                                                                                                              + rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_N )])
                                                                                                                + 0.5 * 0.5 * ( rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C )]
                                                                                                                                + rhsP1[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_E )]));
618
619
    }

620
    dst[vertexdof::macroedge::indexFromVertex( level, i, sD::VERTEX_C )] = ValueType( tmp );
621
622
623
  }
}

624
template< typename ValueType >
625
inline void printFunctionMemory( const uint_t & level, const Edge& edge, const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &dstId){
626
  ValueType* edgeMemory = edge.getData(dstId)->getPointer( level );
Dominik Thoennes's avatar
Dominik Thoennes committed
627
628
629
  using namespace std;
  cout << setfill('=') << setw(100) << "" << endl;
  cout << edge << std::left << setprecision(1) << fixed << setfill(' ') << endl;
630
  uint_t rowsize = levelinfo::num_microvertices_per_edge( level );
Dominik Thoennes's avatar
Dominik Thoennes committed
631
632
633

  if(edge.getNumNeighborFaces() == 2) {
    for (uint_t i = 0; i < (rowsize - 1); ++i) {
634
      cout << setw(5) << edgeMemory[hyteg::vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_N )] << "|";
Dominik Thoennes's avatar
Dominik Thoennes committed
635
636
637
638
    }
    cout << endl;
  }
  for(uint_t i = 0; i < rowsize; ++i){
639
    cout << setw(5) << edgeMemory[hyteg::vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )] << "|";
Dominik Thoennes's avatar
Dominik Thoennes committed
640
641
642
  }
  cout << endl << "     |";
  for(uint_t i = 1; i < rowsize; ++i){
643
    cout << setw(5) << edgeMemory[hyteg::vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_S )] << "|";
Dominik Thoennes's avatar
Dominik Thoennes committed
644
645
646
647
648
649
  }
  cout << endl << setfill('=') << setw(100) << "" << endl << setfill(' ');

}


650
template< typename ValueType >
Dominik Thoennes's avatar
Dominik Thoennes committed
651
inline ValueType getMaxValue( const uint_t & level, Edge &edge, const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &srcId ) {
652
653
654
655

  uint_t rowsize = levelinfo::num_microvertices_per_edge( level );

  auto src = edge.getData( srcId )->getPointer( level );
Dominik Thoennes's avatar
Dominik Thoennes committed
656
  auto localMax = -std::numeric_limits< ValueType >::max();
657
658
659
660
661
662
663
664
665
666

  for( size_t i = 1; i < rowsize - 1; ++i ) {
    localMax = std::max( localMax, src[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C ) ] );
  }

  return localMax;
}


template< typename ValueType >
Dominik Thoennes's avatar
Dominik Thoennes committed
667
inline ValueType getMaxMagnitude( const uint_t & level, Edge &edge, const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &srcId ) {
668
669
670
671

  uint_t rowsize = levelinfo::num_microvertices_per_edge( level );

  auto src = edge.getData( srcId )->getPointer( level );
Dominik Thoennes's avatar
Dominik Thoennes committed
672
  auto localMax = ValueType(0.0);
673
674
675
676
677
678
679
680
681

  for( size_t i = 1; i < rowsize - 1; ++i ) {
    localMax = std::max( localMax, std::abs( src[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C ) ] ));
  }

  return localMax;
}

template< typename ValueType >
Dominik Thoennes's avatar
Dominik Thoennes committed
682
inline ValueType getMinValue( const uint_t & level, Edge &edge, const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &srcId ) {
683
684
685
686

  uint_t rowsize = levelinfo::num_microvertices_per_edge( level );

  auto src = edge.getData( srcId )->getPointer( level );
Dominik Thoennes's avatar
Dominik Thoennes committed
687
  auto localMin = std::numeric_limits< ValueType >::max();
688
689
690
691
692
693
694
695
696

  for( size_t i = 1; i < rowsize - 1; ++i ) {
    localMin = std::min( localMin, src[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C ) ] );
  }

  return localMin;
}


697

Dominik Thoennes's avatar
Dominik Thoennes committed
698
#ifdef HYTEG_BUILD_WITH_PETSC
699
700
701
702
703
704
705
inline void saveOperator( const uint_t&                                           level,
                          Edge&                                                   edge,
                          const PrimitiveStorage&                                 storage,
                          const PrimitiveDataID< StencilMemory< real_t >, Edge >& operatorId,
                          const PrimitiveDataID< FunctionMemory< idx_t >, Edge >& srcId,
                          const PrimitiveDataID< FunctionMemory< idx_t >, Edge >& dstId,
                          const std::shared_ptr< SparseMatrixProxy >&             mat )
706
{
707
   size_t rowsize = levelinfo::num_microvertices_per_edge( level );
708

709
710
711
   auto opr_data = edge.getData( operatorId )->getPointer( level );
   auto src      = edge.getData( srcId )->getPointer( level );
   auto dst      = edge.getData( dstId )->getPointer( level );
712

713
714
   for ( uint_t i = 1; i < rowsize - 1; ++i )
   {
715
716
      idx_t dstint = dst[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )];
      idx_t srcint = src[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )];
717
718
      mat->addValue(
          uint_c( dstint ), uint_c( srcint ), opr_data[vertexdof::stencilIndexFromVertex( stencilDirection::VERTEX_C )] );
719

720
721
722
723
724
      for ( const auto& neighbor : vertexdof::macroedge::neighborsOnEdgeFromVertexDoF )
      {
         srcint = src[vertexdof::macroedge::indexFromVertex( level, i, neighbor )];
         mat->addValue( uint_c( dstint ), uint_c( srcint ), opr_data[vertexdof::stencilIndexFromVertex( neighbor )] );
      }
725

726
727
728
729
730
731
732
733
734
735
736
737
      for ( uint_t neighborFace = 0; neighborFace < edge.getNumNeighborFaces(); neighborFace++ )
      {
         srcint = src[vertexdof::macroedge::indexFromVertexOnNeighborFace( level, i, neighborFace, stencilDirection::VERTEX_W )];
         mat->addValue( uint_c( dstint ),
                        uint_c( srcint ),
                        opr_data[vertexdof::macroedge::stencilIndexOnNeighborFace( stencilDirection::VERTEX_W, neighborFace )] );

         srcint = src[vertexdof::macroedge::indexFromVertexOnNeighborFace( level, i, neighborFace, stencilDirection::VERTEX_E )];
         mat->addValue( uint_c( dstint ),
                        uint_c( srcint ),
                        opr_data[vertexdof::macroedge::stencilIndexOnNeighborFace( stencilDirection::VERTEX_E, neighborFace )] );
      }
738

739
      for ( uint_t neighborCellID = 0; neighborCellID < edge.getNumNeighborCells(); neighborCellID++ )
740
      {
741
742
743
744
745
746
747
748
749
750
751
752
753
754
         const auto& neighborCell              = *( storage.getCell( edge.neighborCells().at( neighborCellID ) ) );
         const auto  localEdgeIDOnNeighborCell = neighborCell.getLocalEdgeID( edge.getID() );
         if ( localEdgeIDOnNeighborCell == 1 || localEdgeIDOnNeighborCell == 4 )
         {
            // Since the functions we access here carry the petsc vector indices, we cannot simply also loop over
            // ghost layer DoFs that do not exist. In the apply kernel this is okay, as we only add zeros in that case.
            // Therefore we check if there are inner vertices - this only applies for macro-edge IDs 1 and 4.
            srcint =
                src[vertexdof::macroedge::indexFromVertexOnNeighborCell( level, i, neighborCellID, edge.getNumNeighborFaces() )];
            mat->addValue(
                uint_c( dstint ),
                uint_c( srcint ),
                opr_data[vertexdof::macroedge::stencilIndexOnNeighborCell( neighborCellID, edge.getNumNeighborFaces() )] );
         }
755
      }
756
   }
757
758
}

759
760
761
762
inline void saveIdentityOperator( const uint_t&                                           level,
                                  Edge&                                                   edge,
                                  const PrimitiveDataID< FunctionMemory< idx_t >, Edge >& dstId,
                                  const std::shared_ptr< SparseMatrixProxy >&             mat )
763
764
765
766
767
768
769
{
   size_t rowsize = levelinfo::num_microvertices_per_edge( level );

   auto dst = edge.getData( dstId )->getPointer( level );

   for ( uint_t i = 1; i < rowsize - 1; ++i )
   {
770
      idx_t dstint = dst[vertexdof::macroedge::indexFromVertex( level, i, stencilDirection::VERTEX_C )];
771
772
773
774
      mat->addValue( uint_c( dstint ), uint_c( dstint ), 1.0 );
   }
}

775
776
777
778
template < typename ValueType >
inline void createVectorFromFunction( const uint_t&                                               level,
                                      Edge&                                                       edge,
                                      const PrimitiveDataID< FunctionMemory< ValueType >, Edge >& srcId,
779
                                      const PrimitiveDataID< FunctionMemory< idx_t >, Edge >&     numeratorId,
780
781
                                      const std::shared_ptr< VectorProxy >&                       vec )
{
782
   idx_t rowsize = (idx_t) levelinfo::num_microvertices_per_edge( level );
783

784
785
   auto src       = edge.getData( srcId )->getPointer( level );
   auto numerator = edge.getData( numeratorId )->getPointer( level );
786

787
788
789
790
   for ( uint_t i = 1; i < uint_c( rowsize - 1 ); i++ )
   {
      vec->setValue( numerator[i], src[i] );
   }
791
}
792

793
794
795
796
template < typename ValueType >
inline void createFunctionFromVector( const uint_t&                                               level,
                                      Edge&                                                       edge,
                                      const PrimitiveDataID< FunctionMemory< ValueType >, Edge >& srcId,
797
                                      const PrimitiveDataID< FunctionMemory< idx_t >, Edge >&     numeratorId,
798
799
                                      const std::shared_ptr< VectorProxy >&                       vec )
{
800
   idx_t rowsize = (idx_t) levelinfo::num_microvertices_per_edge( level );
801

802
   auto numerator = edge.getData( numeratorId )->getPointer( level );
803

804
805
806
807
   for ( uint_t i = 1; i < uint_c( rowsize - 1 ); i++ )
   {
      edge.getData( srcId )->getPointer( level )[i] = vec->getValue( numerator[i] );
   }
808
}
809

810
811
inline void applyDirichletBC( const uint_t&                                           level,
                              Edge&                                                   edge,
812
                              std::vector< idx_t >&                                   mat,
813
814
815
                              const PrimitiveDataID< FunctionMemory< idx_t >, Edge >& numeratorId )
{
   size_t rowsize = levelinfo::num_microvertices_per_edge( level );
816

817
818
   for ( uint_t i = 1; i < rowsize - 1; i++ )
   {
819
      mat.push_back( edge.getData( numeratorId )->getPointer( level )[i] );
820
   }
821
}
Daniel Drzisga's avatar
Daniel Drzisga committed
822
#endif
823

824
825
} // namespace macroedge
} // namespace vertexdof
826
} // namespace hyteg