EdgeDoFPetsc.hpp 32.4 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, 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/>.
 */
Daniel Drzisga's avatar
Daniel Drzisga committed
20
21
22
23
24

#pragma once

#include "core/DataTypes.h"

25
#include "hyteg/edgedofspace/EdgeDoFFunction.hpp"
Dominik Thoennes's avatar
Dominik Thoennes committed
26
27
28
#include "hyteg/edgedofspace/EdgeDoFMacroEdge.hpp"
#include "hyteg/edgedofspace/EdgeDoFMacroFace.hpp"
#include "hyteg/petsc/PETScWrapper.hpp"
29
#include "hyteg/sparseassembly/SparseMatrixProxy.hpp"
30
#include "hyteg/sparseassembly/VectorProxy.hpp"
Daniel Drzisga's avatar
Daniel Drzisga committed
31

32
namespace hyteg {
33
namespace edgedof {
Daniel Drzisga's avatar
Daniel Drzisga committed
34
35
36
37

using walberla::real_t;
using walberla::uint_t;

Dominik Thoennes's avatar
Dominik Thoennes committed
38
#ifdef HYTEG_BUILD_WITH_PETSC
Daniel Drzisga's avatar
Daniel Drzisga committed
39

40
41
42
43
44
45
inline void saveEdgeOperator( const uint_t&                                           Level,
                              const Edge&                                             edge,
                              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 )
Daniel Drzisga's avatar
Daniel Drzisga committed
46
{
47
48
49
50
   using namespace hyteg::edgedof::macroedge;
   size_t rowsize = levelinfo::num_microedges_per_edge( Level );

   real_t*   opr_data = edge.getData( operatorId )->getPointer( Level );
51
52
   idx_t*    src      = edge.getData( srcId )->getPointer( Level );
   idx_t*    dst      = edge.getData( dstId )->getPointer( Level );
53

54
55
   idx_t srcInt;
   idx_t dstInt;
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

   for ( uint_t i = 0; i < rowsize; ++i )
   {
      dstInt = dst[indexFromHorizontalEdge( Level, i, stencilDirection::EDGE_HO_C )];

      for ( uint_t k = 0; k < neighborsOnEdgeFromHorizontalEdge.size(); ++k )
      {
         srcInt = src[indexFromHorizontalEdge( Level, i, neighborsOnEdgeFromHorizontalEdge[k] )];
         mat->addValue( uint_c( dstInt ),
                        uint_c( srcInt ),
                        opr_data[hyteg::edgedof::stencilIndexFromHorizontalEdge( neighborsOnEdgeFromHorizontalEdge[k] )] );
      }
      for ( uint_t k = 0; k < neighborsOnSouthFaceFromHorizontalEdge.size(); ++k )
      {
         srcInt = src[indexFromHorizontalEdge( Level, i, neighborsOnSouthFaceFromHorizontalEdge[k] )];
         mat->addValue( uint_c( dstInt ),
                        uint_c( srcInt ),
                        opr_data[hyteg::edgedof::stencilIndexFromHorizontalEdge( neighborsOnSouthFaceFromHorizontalEdge[k] )] );
Daniel Drzisga's avatar
Daniel Drzisga committed
74
      }
75
76
77
78
79
80
81
82
83
84
85
86
      if ( edge.getNumNeighborFaces() == 2 )
      {
         for ( uint_t k = 0; k < neighborsOnNorthFaceFromHorizontalEdge.size(); ++k )
         {
            srcInt = src[indexFromHorizontalEdge( Level, i, neighborsOnNorthFaceFromHorizontalEdge[k] )];
            mat->addValue(
                uint_c( dstInt ),
                uint_c( srcInt ),
                opr_data[hyteg::edgedof::stencilIndexFromHorizontalEdge( neighborsOnNorthFaceFromHorizontalEdge[k] )] );
         }
      }
   }
Daniel Drzisga's avatar
Daniel Drzisga committed
87
88
}

89
90
91
92
inline void saveEdgeIdentityOperator( const uint_t&                                           Level,
                                      const Edge&                                             edge,
                                      const PrimitiveDataID< FunctionMemory< idx_t >, Edge >& dstId,
                                      const std::shared_ptr< SparseMatrixProxy >&             mat )
93
94
95
96
{
   using namespace hyteg::edgedof::macroedge;
   size_t rowsize = levelinfo::num_microedges_per_edge( Level );

97
98
   idx_t* dst = edge.getData( dstId )->getPointer( Level );
   idx_t  dstInt;
99
100
101
102
103
104
105

   for ( uint_t i = 0; i < rowsize; ++i )
   {
      dstInt = dst[indexFromHorizontalEdge( Level, i, stencilDirection::EDGE_HO_C )];
      mat->addValue( uint_c( dstInt ), uint_c( dstInt ), 1.0 );
   }
}
Daniel Drzisga's avatar
Daniel Drzisga committed
106

107
108
109
110
inline void saveEdgeOperator3D( const uint_t&                                                              level,
                                const Edge&                                                                edge,
                                const PrimitiveStorage&                                                    storage,
                                const PrimitiveDataID< LevelWiseMemory< macroedge::StencilMap_T >, Edge >& operatorId,
111
112
                                const PrimitiveDataID< FunctionMemory< idx_t >, Edge >&                    srcId,
                                const PrimitiveDataID< FunctionMemory< idx_t >, Edge >&                    dstId,
113
                                const std::shared_ptr< SparseMatrixProxy >&                                mat )
Nils Kohl's avatar
Nils Kohl committed
114
{
115
116
117
   auto opr_data = edge.getData( operatorId )->getData( level );
   auto src      = edge.getData( srcId )->getPointer( level );
   auto dst      = edge.getData( dstId )->getPointer( level );
Nils Kohl's avatar
Nils Kohl committed
118

119
120
121
   for ( const auto& centerIndexOnEdge : hyteg::edgedof::macroedge::Iterator( level, 0 ) )
   {
      const EdgeDoFOrientation edgeCenterOrientation = EdgeDoFOrientation::X;
Nils Kohl's avatar
Nils Kohl committed
122

123
      const auto dstInt = dst[edgedof::macroedge::index( level, centerIndexOnEdge.x() )];
Nils Kohl's avatar
Nils Kohl committed
124

125
      for ( uint_t neighborCellID = 0; neighborCellID < edge.getNumNeighborCells(); neighborCellID++ )
Nils Kohl's avatar
Nils Kohl committed
126
      {
127
128
         const Cell& neighborCell    = *( storage.getCell( edge.neighborCells().at( neighborCellID ) ) );
         auto        cellLocalEdgeID = neighborCell.getLocalEdgeID( edge.getID() );
Nils Kohl's avatar
Nils Kohl committed
129

130
131
132
         const auto basisInCell = algorithms::getMissingIntegersAscending< 2, 4 >(
             {neighborCell.getEdgeLocalVertexToCellLocalVertexMaps().at( cellLocalEdgeID ).at( 0 ),
              neighborCell.getEdgeLocalVertexToCellLocalVertexMaps().at( cellLocalEdgeID ).at( 1 )} );
Nils Kohl's avatar
Nils Kohl committed
133

134
135
136
137
         const auto centerIndexInCell = indexing::basisConversion(
             centerIndexOnEdge, basisInCell, {0, 1, 2, 3}, levelinfo::num_microedges_per_edge( level ) );
         const auto cellCenterOrientation = edgedof::convertEdgeDoFOrientationFaceToCell(
             edgeCenterOrientation, basisInCell.at( 0 ), basisInCell.at( 1 ), basisInCell.at( 2 ) );
Nils Kohl's avatar
Nils Kohl committed
138

139
140
141
         for ( const auto& leafOrientationInCell : edgedof::allEdgeDoFOrientations )
         {
            for ( const auto& stencilIt : opr_data[neighborCellID][cellCenterOrientation][leafOrientationInCell] )
Nils Kohl's avatar
Nils Kohl committed
142
            {
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
               const auto stencilOffset = stencilIt.first;
               const auto stencilWeight = stencilIt.second;

               const auto leafOrientationOnEdge = edgedof::convertEdgeDoFOrientationCellToFace(
                   leafOrientationInCell, basisInCell.at( 0 ), basisInCell.at( 1 ), basisInCell.at( 2 ) );
               const auto leafIndexInCell = centerIndexInCell + stencilOffset;

               const auto leafIndexOnEdge = indexing::basisConversion(
                   leafIndexInCell, {0, 1, 2, 3}, basisInCell, levelinfo::num_microedges_per_edge( level ) );

               const auto onCellFacesSet = edgedof::macrocell::isOnCellFaces( level, leafIndexInCell, leafOrientationInCell );
               const auto onCellFacesSetOnEdge =
                   edgedof::macrocell::isOnCellFaces( level, leafIndexOnEdge, leafOrientationOnEdge );

               WALBERLA_ASSERT_EQUAL( onCellFacesSet.size(), onCellFacesSetOnEdge.size() );

               uint_t leafArrayIndexOnEdge = std::numeric_limits< uint_t >::max();

               const auto cellLocalIDsOfNeighborFaces =
                   indexing::cellLocalEdgeIDsToCellLocalNeighborFaceIDs.at( cellLocalEdgeID );
               std::vector< uint_t > cellLocalIDsOfNeighborFacesWithLeafOnThem;
               std::set_intersection( cellLocalIDsOfNeighborFaces.begin(),
                                      cellLocalIDsOfNeighborFaces.end(),
                                      onCellFacesSet.begin(),
                                      onCellFacesSet.end(),
                                      std::back_inserter( cellLocalIDsOfNeighborFacesWithLeafOnThem ) );

               if ( cellLocalIDsOfNeighborFacesWithLeafOnThem.size() == 0 )
               {
                  // leaf in macro-cell
                  leafArrayIndexOnEdge = edgedof::macroedge::indexOnNeighborCell(
                      level, leafIndexOnEdge.x(), neighborCellID, edge.getNumNeighborFaces(), leafOrientationOnEdge );
               }
               else if ( cellLocalIDsOfNeighborFacesWithLeafOnThem.size() == 1 )
               {
                  // leaf on macro-face
                  WALBERLA_ASSERT( !edgedof::macrocell::isInnerEdgeDoF( level, leafIndexInCell, leafOrientationInCell ) );

                  const auto cellLocalFaceID = *cellLocalIDsOfNeighborFacesWithLeafOnThem.begin();
                  const auto facePrimitiveID = neighborCell.neighborFaces().at( cellLocalFaceID );
                  WALBERLA_ASSERT( std::find( edge.neighborFaces().begin(), edge.neighborFaces().end(), facePrimitiveID ) !=
                                   edge.neighborFaces().end() );

                  // The leaf orientation on the edge must be X, Y or XY since it is located on a neighboring face.
                  // Therefore we need to know the three spanning vertex IDs and convert the leaf orientation again.
                  const auto spanningCellLocalVertices = indexing::cellLocalFaceIDsToSpanningVertexIDs.at( cellLocalFaceID );
                  std::array< uint_t, 4 > faceBasisInCell;
                  if ( spanningCellLocalVertices.count( basisInCell.at( 2 ) ) == 1 )
                  {
                     faceBasisInCell = basisInCell;
                  }
                  else
                  {
                     WALBERLA_ASSERT( spanningCellLocalVertices.count( basisInCell.at( 3 ) ) == 1 );
                     faceBasisInCell    = basisInCell;
                     faceBasisInCell[2] = basisInCell.at( 3 );
                     faceBasisInCell[3] = basisInCell.at( 2 );
                  }

                  const auto leafIndexOnEdgeGhostLayer = indexing::basisConversion(
                      leafIndexInCell, {0, 1, 2, 3}, faceBasisInCell, levelinfo::num_microedges_per_edge( level ) );
                  const auto leafOrientationOnEdgeGhostLayer = edgedof::convertEdgeDoFOrientationCellToFace(
                      leafOrientationInCell, faceBasisInCell.at( 0 ), faceBasisInCell.at( 1 ), faceBasisInCell.at( 2 ) );

                  const auto localFaceIDOnEdge = edge.face_index( facePrimitiveID );
                  leafArrayIndexOnEdge         = edgedof::macroedge::indexOnNeighborFace(
                      level, leafIndexOnEdgeGhostLayer.x(), localFaceIDOnEdge, leafOrientationOnEdgeGhostLayer );
               }
               else
               {
                  // leaf on macro-edge
                  WALBERLA_ASSERT_EQUAL( cellLocalIDsOfNeighborFacesWithLeafOnThem.size(), 2 );
                  WALBERLA_ASSERT( !edgedof::macrocell::isInnerEdgeDoF( level, leafIndexInCell, leafOrientationInCell ) );
                  WALBERLA_ASSERT_EQUAL( leafOrientationOnEdge, EdgeDoFOrientation::X );
                  leafArrayIndexOnEdge = edgedof::macroedge::index( level, leafIndexOnEdge.x() );
               }

               const auto srcInt = src[leafArrayIndexOnEdge];
               mat->addValue( uint_c( dstInt ), uint_c( srcInt ), stencilWeight );
Nils Kohl's avatar
Nils Kohl committed
222
            }
223
         }
Nils Kohl's avatar
Nils Kohl committed
224
      }
225
   }
Nils Kohl's avatar
Nils Kohl committed
226
227
}

228
229
230
231
232
233
inline void saveFaceOperator( const uint_t&                                           Level,
                              const Face&                                             face,
                              const PrimitiveDataID< StencilMemory< real_t >, Face >& operatorId,
                              const PrimitiveDataID< FunctionMemory< idx_t >, Face >& srcId,
                              const PrimitiveDataID< FunctionMemory< idx_t >, Face >& dstId,
                              const std::shared_ptr< SparseMatrixProxy >&             mat )
Daniel Drzisga's avatar
Daniel Drzisga committed
234
{
235
   real_t*   opr_data = face.getData( operatorId )->getPointer( Level );
236
237
   idx_t*    src      = face.getData( srcId )->getPointer( Level );
   idx_t*    dst      = face.getData( dstId )->getPointer( Level );
238

239
240
   idx_t srcInt;
   idx_t dstInt;
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255

   using namespace edgedof::macroface;

   for ( const auto& it : hyteg::edgedof::macroface::Iterator( Level, 0 ) )
   {
      if ( it.row() != 0 )
      {
         dstInt = dst[indexFromHorizontalEdge( Level, it.col(), it.row(), stencilDirection::EDGE_HO_C )];
         for ( uint_t k = 0; k < neighborsFromHorizontalEdge.size(); ++k )
         {
            srcInt = src[indexFromHorizontalEdge( Level, it.col(), it.row(), neighborsFromHorizontalEdge[k] )];
            mat->addValue( uint_c( dstInt ),
                           uint_c( srcInt ),
                           opr_data[edgedof::stencilIndexFromHorizontalEdge( neighborsFromHorizontalEdge[k] )] );
         }
Daniel Drzisga's avatar
Daniel Drzisga committed
256
      }
257
258
259
260
261
262
263
264
265
266
      if ( it.col() + it.row() != ( hyteg::levelinfo::num_microedges_per_edge( Level ) - 1 ) )
      {
         dstInt = dst[indexFromDiagonalEdge( Level, it.col(), it.row(), stencilDirection::EDGE_DI_C )];
         for ( uint_t k = 0; k < neighborsFromDiagonalEdge.size(); ++k )
         {
            srcInt = src[indexFromDiagonalEdge( Level, it.col(), it.row(), neighborsFromDiagonalEdge[k] )];
            mat->addValue( uint_c( dstInt ),
                           uint_c( srcInt ),
                           opr_data[edgedof::stencilIndexFromDiagonalEdge( neighborsFromDiagonalEdge[k] )] );
         }
Daniel Drzisga's avatar
Daniel Drzisga committed
267
      }
268
269
270
271
272
273
274
275
276
277
      if ( it.col() != 0 )
      {
         dstInt = dst[indexFromVerticalEdge( Level, it.col(), it.row(), stencilDirection::EDGE_VE_C )];
         for ( uint_t k = 0; k < neighborsFromVerticalEdge.size(); ++k )
         {
            srcInt = src[indexFromVerticalEdge( Level, it.col(), it.row(), neighborsFromVerticalEdge[k] )];
            mat->addValue( uint_c( dstInt ),
                           uint_c( srcInt ),
                           opr_data[edgedof::stencilIndexFromVerticalEdge( neighborsFromVerticalEdge[k] )] );
         }
Daniel Drzisga's avatar
Daniel Drzisga committed
278
      }
279
   }
Daniel Drzisga's avatar
Daniel Drzisga committed
280
281
}

282
283
284
285
inline void saveFaceIdentityOperator( const uint_t&                                           Level,
                                      const Face&                                             face,
                                      const PrimitiveDataID< FunctionMemory< idx_t >, Face >& dstId,
                                      const std::shared_ptr< SparseMatrixProxy >&             mat )
286
{
287
   idx_t* dst = face.getData( dstId )->getPointer( Level );
288

289
   idx_t dstInt;
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312

   using namespace edgedof::macroface;

   for ( const auto& it : hyteg::edgedof::macroface::Iterator( Level, 0 ) )
   {
      if ( it.row() != 0 )
      {
         dstInt = dst[indexFromHorizontalEdge( Level, it.col(), it.row(), stencilDirection::EDGE_HO_C )];
         mat->addValue( uint_c( dstInt ), uint_c( dstInt ), 1.0 );
      }
      if ( it.col() + it.row() != ( hyteg::levelinfo::num_microedges_per_edge( Level ) - 1 ) )
      {
         dstInt = dst[indexFromDiagonalEdge( Level, it.col(), it.row(), stencilDirection::EDGE_DI_C )];
         mat->addValue( uint_c( dstInt ), uint_c( dstInt ), 1.0 );
      }
      if ( it.col() != 0 )
      {
         dstInt = dst[indexFromVerticalEdge( Level, it.col(), it.row(), stencilDirection::EDGE_VE_C )];
         mat->addValue( uint_c( dstInt ), uint_c( dstInt ), 1.0 );
      }
   }
}

313
314
315
316
inline void saveFaceOperator3D( const uint_t&                                                              level,
                                const Face&                                                                face,
                                const PrimitiveStorage&                                                    storage,
                                const PrimitiveDataID< LevelWiseMemory< macroface::StencilMap_T >, Face >& operatorId,
317
318
                                const PrimitiveDataID< FunctionMemory< idx_t >, Face >&                    srcId,
                                const PrimitiveDataID< FunctionMemory< idx_t >, Face >&                    dstId,
319
                                const std::shared_ptr< SparseMatrixProxy >&                                mat )
320
{
321
322
323
   auto opr_data = face.getData( operatorId )->getData( level );
   auto src      = face.getData( srcId )->getPointer( level );
   auto dst      = face.getData( dstId )->getPointer( level );
324

325
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
   for ( const auto& centerIndexInFace : hyteg::edgedof::macroface::Iterator( level, 0 ) )
   {
      for ( const auto& faceCenterOrientation : edgedof::faceLocalEdgeDoFOrientations )
      {
         if ( faceCenterOrientation == edgedof::EdgeDoFOrientation::X &&
              edgedof::macroface::isHorizontalEdgeOnBoundary( level, centerIndexInFace ) )
            continue;
         if ( faceCenterOrientation == edgedof::EdgeDoFOrientation::Y &&
              edgedof::macroface::isVerticalEdgeOnBoundary( level, centerIndexInFace ) )
            continue;
         if ( faceCenterOrientation == edgedof::EdgeDoFOrientation::XY &&
              edgedof::macroface::isDiagonalEdgeOnBoundary( level, centerIndexInFace ) )
            continue;

         const auto dstIdx =
             edgedof::macroface::index( level, centerIndexInFace.x(), centerIndexInFace.y(), faceCenterOrientation );
         const auto dstInt = dst[dstIdx];

         for ( uint_t neighborCellID = 0; neighborCellID < face.getNumNeighborCells(); neighborCellID++ )
         {
            const Cell&  neighborCell = *( storage.getCell( face.neighborCells().at( neighborCellID ) ) );
            const uint_t localFaceID  = neighborCell.getLocalFaceID( face.getID() );

            const auto centerIndexInCell =
                macroface::getIndexInNeighboringMacroCell( centerIndexInFace, face, neighborCellID, storage, level );
            const auto cellCenterOrientation =
                macroface::getOrientattionInNeighboringMacroCell( faceCenterOrientation, face, neighborCellID, storage );

            for ( const auto& leafOrientation : edgedof::allEdgeDoFOrientations )
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
384
385
386
387
               for ( const auto& stencilIt : opr_data[neighborCellID][cellCenterOrientation][leafOrientation] )
               {
                  const auto stencilOffset = stencilIt.first;
                  const auto stencilWeight = stencilIt.second;

                  const auto leafOrientationInFace =
                      macrocell::getOrientattionInNeighboringMacroFace( leafOrientation, neighborCell, localFaceID, storage );

                  const auto leafIndexInCell = centerIndexInCell + stencilOffset;
                  const auto leafIndexInFace =
                      leafOrientation == edgedof::EdgeDoFOrientation::XYZ ?
                          macrocell::getIndexInNeighboringMacroFaceXYZ(
                              leafIndexInCell, neighborCell, localFaceID, storage, level ) :
                          macrocell::getIndexInNeighboringMacroFace( leafIndexInCell, neighborCell, localFaceID, storage, level );

                  WALBERLA_ASSERT_LESS_EQUAL( leafIndexInFace.z(), 1 );

                  uint_t leafArrayIndexInFace;
                  if ( algorithms::contains( edgedof::faceLocalEdgeDoFOrientations, leafOrientationInFace ) &&
                       leafIndexInFace.z() == 0 )
                  {
                     leafArrayIndexInFace =
                         edgedof::macroface::index( level, leafIndexInFace.x(), leafIndexInFace.y(), leafOrientationInFace );
                  }
                  else
                  {
                     leafArrayIndexInFace = edgedof::macroface::index(
                         level, leafIndexInFace.x(), leafIndexInFace.y(), leafOrientationInFace, neighborCellID );
                  }

                  const auto srcInt = src[leafArrayIndexInFace];
                  mat->addValue( uint_c( dstInt ), uint_c( srcInt ), stencilWeight );
               }
388
            }
389
         }
390
      }
391
   }
392
393
}

394
395
396
inline void saveCellOperator( const uint_t&                                                              Level,
                              const Cell&                                                                cell,
                              const PrimitiveDataID< LevelWiseMemory< macrocell::StencilMap_T >, Cell >& operatorId,
397
398
                              const PrimitiveDataID< FunctionMemory< idx_t >, Cell >&                    srcId,
                              const PrimitiveDataID< FunctionMemory< idx_t >, Cell >&                    dstId,
399
                              const std::shared_ptr< SparseMatrixProxy >&                                mat )
400
{
401
402
403
404
405
406
407
   auto srcData  = cell.getData( srcId )->getPointer( Level );
   auto dstData  = cell.getData( dstId )->getPointer( Level );
   auto opr_data = cell.getData( operatorId )->getData( Level );

   for ( const auto& it : edgedof::macrocell::Iterator( Level, 0 ) )
   {
      std::vector< edgedof::EdgeDoFOrientation > innerOrientations;
408

409
410
411
412
413
414
415
416
417
418
419
420
421
422
      if ( macrocell::isInnerXEdgeDoF( Level, it ) )
         innerOrientations.push_back( edgedof::EdgeDoFOrientation::X );
      if ( macrocell::isInnerYEdgeDoF( Level, it ) )
         innerOrientations.push_back( edgedof::EdgeDoFOrientation::Y );
      if ( macrocell::isInnerZEdgeDoF( Level, it ) )
         innerOrientations.push_back( edgedof::EdgeDoFOrientation::Z );
      if ( macrocell::isInnerXYEdgeDoF( Level, it ) )
         innerOrientations.push_back( edgedof::EdgeDoFOrientation::XY );
      if ( macrocell::isInnerXZEdgeDoF( Level, it ) )
         innerOrientations.push_back( edgedof::EdgeDoFOrientation::XZ );
      if ( macrocell::isInnerYZEdgeDoF( Level, it ) )
         innerOrientations.push_back( edgedof::EdgeDoFOrientation::YZ );

      for ( const auto& centerOrientation : innerOrientations )
423
      {
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
         const auto dstArrayIdx = edgedof::macrocell::index( Level, it.x(), it.y(), it.z(), centerOrientation );
         const auto dstInt      = dstData[dstArrayIdx];

         for ( const auto& leafOrientation : edgedof::allEdgeDoFOrientations )
         {
            const auto edgeDoFNeighbors =
                P2Elements::P2Elements3D::getAllEdgeDoFNeighborsFromEdgeDoFInMacroCell( centerOrientation, leafOrientation );
            for ( const auto& neighbor : edgeDoFNeighbors )
            {
               const auto   srcIdx      = it + neighbor;
               const auto   srcArrayIdx = edgedof::macrocell::index( Level, srcIdx.x(), srcIdx.y(), srcIdx.z(), leafOrientation );
               const auto   srcInt      = srcData[srcArrayIdx];
               const real_t stencilWeight = opr_data[centerOrientation][leafOrientation][neighbor];
               mat->addValue( uint_c( dstInt ), uint_c( srcInt ), stencilWeight );
            }
         }
440
      }
441
   }
442

443
444
445
   for ( const auto& it : edgedof::macrocell::IteratorXYZ( Level, 0 ) )
   {
      const auto centerOrientation = edgedof::EdgeDoFOrientation::XYZ;
446

447
448
      const auto dstArrayIdx = edgedof::macrocell::index( Level, it.x(), it.y(), it.z(), centerOrientation );
      const auto dstInt      = dstData[dstArrayIdx];
449

450
      for ( const auto& leafOrientation : edgedof::allEdgeDoFOrientations )
451
      {
452
453
454
455
456
457
458
459
460
461
         const auto edgeDoFNeighbors =
             P2Elements::P2Elements3D::getAllEdgeDoFNeighborsFromEdgeDoFInMacroCell( centerOrientation, leafOrientation );
         for ( const auto& neighbor : edgeDoFNeighbors )
         {
            const auto   srcIdx        = it + neighbor;
            const auto   srcArrayIdx   = edgedof::macrocell::index( Level, srcIdx.x(), srcIdx.y(), srcIdx.z(), leafOrientation );
            const auto   srcInt        = srcData[srcArrayIdx];
            const real_t stencilWeight = opr_data[centerOrientation][leafOrientation][neighbor];
            mat->addValue( uint_c( dstInt ), uint_c( srcInt ), stencilWeight );
         }
462
      }
463
   }
464
465
}

466
467
468
469
inline void saveCellIdentityOperator( const uint_t&                                           Level,
                                      const Cell&                                             cell,
                                      const PrimitiveDataID< FunctionMemory< idx_t >, Cell >& dstId,
                                      const std::shared_ptr< SparseMatrixProxy >&             mat )
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
{
   auto dstData = cell.getData( dstId )->getPointer( Level );

   for ( const auto& it : edgedof::macrocell::Iterator( Level, 0 ) )
   {
      std::vector< edgedof::EdgeDoFOrientation > innerOrientations;

      if ( macrocell::isInnerXEdgeDoF( Level, it ) )
         innerOrientations.push_back( edgedof::EdgeDoFOrientation::X );
      if ( macrocell::isInnerYEdgeDoF( Level, it ) )
         innerOrientations.push_back( edgedof::EdgeDoFOrientation::Y );
      if ( macrocell::isInnerZEdgeDoF( Level, it ) )
         innerOrientations.push_back( edgedof::EdgeDoFOrientation::Z );
      if ( macrocell::isInnerXYEdgeDoF( Level, it ) )
         innerOrientations.push_back( edgedof::EdgeDoFOrientation::XY );
      if ( macrocell::isInnerXZEdgeDoF( Level, it ) )
         innerOrientations.push_back( edgedof::EdgeDoFOrientation::XZ );
      if ( macrocell::isInnerYZEdgeDoF( Level, it ) )
         innerOrientations.push_back( edgedof::EdgeDoFOrientation::YZ );

      for ( const auto& centerOrientation : innerOrientations )
      {
         const auto dstArrayIdx = edgedof::macrocell::index( Level, it.x(), it.y(), it.z(), centerOrientation );
         const auto dstInt      = dstData[dstArrayIdx];
         mat->addValue( uint_c( dstInt ), uint_c( dstInt ), 1.0 );
      }
   }

   for ( const auto& it : edgedof::macrocell::IteratorXYZ( Level, 0 ) )
   {
      const auto centerOrientation = edgedof::EdgeDoFOrientation::XYZ;

      const auto dstArrayIdx = edgedof::macrocell::index( Level, it.x(), it.y(), it.z(), centerOrientation );
      const auto dstInt      = dstData[dstArrayIdx];
      mat->addValue( uint_c( dstInt ), uint_c( dstInt ), 1.0 );
   }
}

508
509
510
511
512
} // namespace edgedof

namespace petsc {

inline void createVectorFromFunction( const EdgeDoFFunction< PetscReal >&   function,
513
                                      const EdgeDoFFunction< idx_t >&       numerator,
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
                                      const std::shared_ptr< VectorProxy >& vec,
                                      uint_t                                level,
                                      DoFType                               flag )
{
   for ( auto& it : function.getStorage()->getEdges() )
   {
      Edge& edge = *it.second;

      const DoFType edgeBC = function.getBoundaryCondition().getBoundaryType( edge.getMeshBoundaryFlag() );
      if ( testFlag( edgeBC, flag ) )
      {
         edgedof::macroedge::createVectorFromFunction< PetscReal >(
             level, edge, function.getEdgeDataID(), numerator.getEdgeDataID(), vec );
      }
   }

   for ( auto& it : function.getStorage()->getFaces() )
   {
      Face& face = *it.second;

      const DoFType faceBC = function.getBoundaryCondition().getBoundaryType( face.getMeshBoundaryFlag() );
      if ( testFlag( faceBC, flag ) )
      {
         edgedof::macroface::createVectorFromFunction< PetscReal >(
             level, face, function.getFaceDataID(), numerator.getFaceDataID(), vec );
      }
   }

   for ( auto& it : function.getStorage()->getCells() )
   {
      Cell& cell = *it.second;

      const DoFType cellBC = function.getBoundaryCondition().getBoundaryType( cell.getMeshBoundaryFlag() );
      if ( testFlag( cellBC, flag ) )
      {
         edgedof::macrocell::createVectorFromFunction< PetscReal >(
             level, cell, function.getCellDataID(), numerator.getCellDataID(), vec );
      }
   }
}

inline void createFunctionFromVector( const EdgeDoFFunction< PetscReal >&   function,
556
                                      const EdgeDoFFunction< idx_t >&       numerator,
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
                                      const std::shared_ptr< VectorProxy >& vec,
                                      uint_t                                level,
                                      DoFType                               flag )
{
   function.startCommunication< Vertex, Edge >( level );
   function.endCommunication< Vertex, Edge >( level );

   for ( auto& it : function.getStorage()->getEdges() )
   {
      Edge& edge = *it.second;

      const DoFType edgeBC = function.getBoundaryCondition().getBoundaryType( edge.getMeshBoundaryFlag() );
      if ( testFlag( edgeBC, flag ) )
      {
         edgedof::macroedge::createFunctionFromVector< PetscReal >(
             level, edge, function.getEdgeDataID(), numerator.getEdgeDataID(), vec );
      }
   }

   function.startCommunication< Edge, Face >( level );
   function.endCommunication< Edge, Face >( level );

   for ( auto& it : function.getStorage()->getFaces() )
   {
      Face& face = *it.second;

      const DoFType faceBC = function.getBoundaryCondition().getBoundaryType( face.getMeshBoundaryFlag() );
      if ( testFlag( faceBC, flag ) )
      {
         edgedof::macroface::createFunctionFromVector< PetscReal >(
             level, face, function.getFaceDataID(), numerator.getFaceDataID(), vec );
      }
   }

   for ( auto& it : function.getStorage()->getCells() )
   {
      Cell& cell = *it.second;

      const DoFType cellBC = function.getBoundaryCondition().getBoundaryType( cell.getMeshBoundaryFlag() );
      if ( testFlag( cellBC, flag ) )
      {
         edgedof::macrocell::createFunctionFromVector< PetscReal >(
             level, cell, function.getCellDataID(), numerator.getCellDataID(), vec );
      }
   }
}

604
inline void applyDirichletBC( const EdgeDoFFunction< idx_t >& numerator, std::vector< idx_t >& mat, uint_t level )
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
{
   for ( auto& it : numerator.getStorage()->getEdges() )
   {
      Edge& edge = *it.second;

      const DoFType edgeBC = numerator.getBoundaryCondition().getBoundaryType( edge.getMeshBoundaryFlag() );
      if ( testFlag( edgeBC, DirichletBoundary ) )
      {
         edgedof::macroedge::applyDirichletBC( level, edge, mat, numerator.getEdgeDataID() );
      }
   }

   for ( auto& it : numerator.getStorage()->getFaces() )
   {
      Face& face = *it.second;

      const DoFType faceBC = numerator.getBoundaryCondition().getBoundaryType( face.getMeshBoundaryFlag() );
      if ( testFlag( faceBC, DirichletBoundary ) )
      {
         edgedof::macroface::applyDirichletBC( level, face, mat, numerator.getFaceDataID() );
      }
   }
}

Dominik Thoennes's avatar
Dominik Thoennes committed
629
template < class OperatorType >
630
inline void createMatrix( const OperatorType&                         opr,
631
632
                          const EdgeDoFFunction< idx_t >&             src,
                          const EdgeDoFFunction< idx_t >&             dst,
633
634
635
                          const std::shared_ptr< SparseMatrixProxy >& mat,
                          size_t                                      level,
                          DoFType                                     flag )
Daniel Drzisga's avatar
Daniel Drzisga committed
636
{
637
   auto storage = src.getStorage();
638

639
640
641
   for ( auto& it : opr.getStorage()->getEdges() )
   {
      Edge& edge = *it.second;
Daniel Drzisga's avatar
Daniel Drzisga committed
642

643
644
      const DoFType edgeBC = dst.getBoundaryCondition().getBoundaryType( edge.getMeshBoundaryFlag() );
      if ( testFlag( edgeBC, flag ) )
Nils Kohl's avatar
Nils Kohl committed
645
      {
646
647
         if ( storage->hasGlobalCells() )
         {
648
649
            edgedof::saveEdgeOperator3D(
                level, edge, *storage, opr.getEdgeStencil3DID(), src.getEdgeDataID(), dst.getEdgeDataID(), mat );
650
651
652
         }
         else
         {
653
            edgedof::saveEdgeOperator( level, edge, opr.getEdgeStencilID(), src.getEdgeDataID(), dst.getEdgeDataID(), mat );
654
         }
Nils Kohl's avatar
Nils Kohl committed
655
      }
656
657
658
659
660
   }

   if ( level >= 1 )
   {
      for ( auto& it : opr.getStorage()->getFaces() )
Nils Kohl's avatar
Nils Kohl committed
661
      {
662
663
664
665
666
667
668
         Face& face = *it.second;

         const DoFType faceBC = dst.getBoundaryCondition().getBoundaryType( face.getMeshBoundaryFlag() );
         if ( testFlag( faceBC, flag ) )
         {
            if ( storage->hasGlobalCells() )
            {
669
               edgedof::saveFaceOperator3D(
670
671
672
673
                   level, face, *storage, opr.getFaceStencil3DID(), src.getFaceDataID(), dst.getFaceDataID(), mat );
            }
            else
            {
674
               edgedof::saveFaceOperator( level, face, opr.getFaceStencilID(), src.getFaceDataID(), dst.getFaceDataID(), mat );
675
676
            }
         }
Nils Kohl's avatar
Nils Kohl committed
677
678
      }

679
680
681
682
683
684
685
      for ( auto& it : opr.getStorage()->getCells() )
      {
         Cell& cell = *it.second;

         const DoFType cellBC = dst.getBoundaryCondition().getBoundaryType( cell.getMeshBoundaryFlag() );
         if ( testFlag( cellBC, flag ) )
         {
686
            edgedof::saveCellOperator( level, cell, opr.getCellStencilID(), src.getCellDataID(), dst.getCellDataID(), mat );
687
688
689
         }
      }
   }
Daniel Drzisga's avatar
Daniel Drzisga committed
690
691
692
693
}

#endif

694
} // namespace edgedof
Daniel Drzisga's avatar
Daniel Drzisga committed
695

696
} // namespace hyteg