40#ifndef DUNE_DIRICHLETCONSTRAINTS_HH
41#define DUNE_DIRICHLETCONSTRAINTS_HH
43#include <dune/fem/function/common/scalarproducts.hh>
44#include <dune/fem/function/localfunction/const.hh>
45#include <dune/fem/operator/common/temporarylocalmatrix.hh>
46#include <dune/fem/common/bindguard.hh>
47#include <dune/fem/common/coordinate.hh>
48#include <dune/fem/common/localcontribution.hh>
49#include <dune/fem/space/common/localinterpolation.hh>
50#include <dune/fem/space/mapper/nonblockmapper.hh>
55class HasLocalFunction;
58template <
class Model,
class DiscreteFunctionSpace >
59class DirichletConstraints
62 enum Operation { set = 0, sub = 1, add = 2 };
63 typedef Model ModelType;
65 typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
66 typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
67 typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;
68 typedef typename DiscreteFunctionSpaceType::HessianRangeType HessianRangeType;
71 typedef typename DiscreteFunctionSpaceType :: GridPartType GridPartType;
73 typedef typename DiscreteFunctionSpaceType :: GridType GridType;
77 static const int localBlockSize = DiscreteFunctionSpaceType :: localBlockSize ;
78 static const int dimRange = DiscreteFunctionSpaceType::FunctionSpaceType::dimRange;
79 static_assert( localBlockSize <= dimRange,
80 "local block size of the space must be less than or equal to the dimension of the range of the function space.");
81 typedef std::array<int,localBlockSize> DirichletBlock;
82 typedef std::vector< DirichletBlock > DirichletBlockVector;
87 typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
89 typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
90 typedef typename DiscreteFunctionSpace::RangeType RangeType;
91 typedef typename DiscreteFunctionSpace::JacobianRangeType JacobianRangeType;
92 typedef typename DiscreteFunctionSpace::HessianRangeType HessianRangeType;
95 const ModelType& impl_;
96 const EntityType& entity_;
101 static const int dimRange = RangeType::dimension;
102 BoundaryWrapper(
const ModelType& impl,
const EntityType& entity,
const int order,
int bndId )
103 : impl_( impl ), entity_(entity), order_(order), bndId_(bndId) {}
104 const EntityType& entity()
const {
return entity_; }
105 const int order ()
const {
return order_; }
106 template <
class Po
int>
107 void evaluate(
const Point& x, RangeType& ret )
const
108 { impl_.dirichlet(bndId_,Dune::Fem::coordinate(x),ret); }
109 template <
class Po
int>
110 void jacobian(
const Point& x, JacobianRangeType& ret )
const
114 DirichletConstraints( ModelType &model,
const DiscreteFunctionSpaceType& space )
119 hasDirichletDofs_( false ),
130 template <
class DiscreteFunctionType >
133 updateDirichletDofs();
136 if( hasDirichletDofs_ )
141 ConstDofIteratorType uIt = u.
dbegin();
142 DofIteratorType wIt = w.
dbegin();
145 const unsigned int blocks = space_.blockMapper().size();
146 for(
unsigned int blockDof = 0; blockDof < blocks ; ++ blockDof )
148 for(
int l = 0; l < localBlockSize ; ++ l, ++ wIt, ++ uIt )
150 if( dirichletBlocks_[ blockDof ][l] )
153 assert( uIt != u.
dend() );
154 assert( wIt != w.
dend() );
168 template <
class DiscreteFunctionType >
171 updateDirichletDofs();
174 if( hasDirichletDofs_ )
178 DofIteratorType wIt = w.
dbegin();
181 const unsigned int blocks = space_.blockMapper().size();
182 for(
unsigned int blockDof = 0; blockDof < blocks ; ++ blockDof )
184 for(
int l = 0; l < localBlockSize ; ++ l, ++ wIt )
186 if( dirichletBlocks_[ blockDof ][l] )
189 assert( wIt != w.
dend() );
203 template <
class DiscreteFunctionType >
206 if (op != Operation::set)
209 return (*
this)(u,w,op);
212 updateDirichletDofs();
214 if( hasDirichletDofs_ )
216 typedef typename DiscreteFunctionSpaceType :: IteratorType IteratorType;
217 typedef typename IteratorType :: Entity EntityType;
219 Dune::Fem::SetSelectedLocalContribution< DiscreteFunctionType > wLocal( w );
220 Dune::Fem::LocalInterpolation< DiscreteFunctionSpaceType > interpolation( space_ );
221 for(
const EntityType &entity : space_ )
226 auto wGuard = Dune::Fem::bindGuard( wLocal, entity );
228 auto iGuard = bindGuard( interpolation, entity );
231 dirichletDofTreatment( interpolation, wLocal );
236 typename = std::enable_if_t< std::is_base_of<Dune::Fem::HasLocalFunction, GridFunctionType>::value > >
237 void operator ()(
const GridFunctionType &u,
240 updateDirichletDofs();
242 if( hasDirichletDofs_ )
244 Dune::Fem::ConstLocalFunction< GridFunctionType > uLocal( u );
245 Dune::Fem::SetSelectedLocalContribution< DiscreteFunctionType > wLocal( w );
246 Dune::Fem::LocalInterpolation< DiscreteFunctionSpaceType > interpolation( space_ );
248 for(
const auto &entity : space_ )
250 auto guard = Dune::Fem::bindGuard( std::tie(uLocal,wLocal), entity );
252 auto iGuard = bindGuard( interpolation, entity );
255 if (op == Operation::sub || op == Operation::add)
257 dirichletDofTreatment( interpolation, uLocal, wLocal, op );
270 template <
class LinearOperator>
271 void applyToOperator( LinearOperator& linearOperator )
const
273 updateDirichletDofs();
275 typedef typename DiscreteFunctionSpaceType :: IteratorType IteratorType;
276 typedef typename IteratorType :: Entity EntityType;
279 if( hasDirichletDofs_ )
281 typedef typename DiscreteFunctionSpaceType :: BlockMapperType BlockMapperType;
284 std::vector<std::size_t> globalBlockDofs;
285 std::vector<std::size_t> globalDofs;
288 std::set<std::size_t> unitRows;
289 std::set<std::size_t> auxRows;
291 const auto& space = space_;
293 const auto &auxiliaryDofs = space.auxiliaryDofs();
295 const IteratorType end = space_.end();
296 for( IteratorType it = space_.begin(); it != end; ++it )
298 const EntityType &entity = *it;
301 const int localBlocks = space.blockMapper().numDofs( entity );
304 globalBlockDofs.resize(localBlocks);
306 space.blockMapper().map( entity, globalBlockDofs );
309 globalDofs.resize(localBlocks * localBlockSize);
310 mapper.map( entity, globalDofs );
322 for(
int localBlockDof = 0 ; localBlockDof < localBlocks; ++ localBlockDof )
324 int global = globalBlockDofs[localBlockDof];
325 std::set<std::size_t>& rows = auxiliaryDofs.contains( global ) ? auxRows : unitRows;
326 for(
int l = 0; l < localBlockSize; ++ l, ++ localDof )
328 if( dirichletBlocks_[global][l] )
331 rows.insert( globalDofs[ localDof ] );
338 linearOperator.setUnitRows( unitRows, auxRows );
342 const DirichletBlockVector &dirichletBlocks()
const
344 updateDirichletDofs();
345 return dirichletBlocks_;
350 template<
class LocalInterpolationType,
class LocalFunctionType >
351 void dirichletDofTreatment(
const LocalInterpolationType& interpolation, LocalFunctionType &wLocal)
const
354 const typename LocalFunctionType::EntityType &entity = wLocal.entity();
357 const int localBlocks = space_.blockMapper().numDofs( entity );
360 std::vector<std::size_t> globalBlockDofs(localBlocks);
361 space_.blockMapper().map( entity, globalBlockDofs );
362 std::vector<typename LocalFunctionType::RangeFieldType> values( localBlocks*localBlockSize );
367 for(
int localBlock = 0 ; localBlock < localBlocks; ++ localBlock )
370 int global = globalBlockDofs[localBlock];
371 for(
int l = 0; l < localBlockSize ; ++ l, ++localDof )
373 if( dirichletBlocks_[ global ][l] )
375 interpolation( BoundaryWrapper(model_, entity, wLocal.order(), dirichletBlocks_[global][l]), values );
377 assert( (
unsigned int)localDof < wLocal.size() );
378 wLocal[ localDof ] = values[ localDof ];
384 template<
class LocalInterpolationType,
class Gr
idLocalFunctionType,
class LocalFunctionType >
385 void dirichletDofTreatment(
const LocalInterpolationType& interpolation,
386 const GridLocalFunctionType &uLocal,
387 LocalFunctionType &wLocal, Operation op )
const
390 const typename LocalFunctionType::EntityType &entity = wLocal.entity();
393 const int localBlocks = space_.blockMapper().numDofs( entity );
395 typedef typename DiscreteFunctionSpaceType::BlockMapperType::GlobalKeyType GlobalKeyType;
398 std::vector< GlobalKeyType > globalBlockDofs(localBlocks);
399 space_.blockMapper().map(entity,globalBlockDofs);
400 std::vector<double> values( localBlocks*localBlockSize );
401 std::vector<double> valuesModel( localBlocks*localBlockSize );
403 interpolation( uLocal, values );
408 for(
int localBlock = 0 ; localBlock < localBlocks; ++ localBlock )
411 int global = globalBlockDofs[localBlock];
412 for(
int l = 0; l < localBlockSize ; ++ l, ++localDof )
414 if( dirichletBlocks_[ global ][l] )
416 if (op == Operation::sub)
418 interpolation(BoundaryWrapper(model_, entity, wLocal.order(), dirichletBlocks_[global][l]), valuesModel);
419 values[ localDof ] -= valuesModel[ localDof ];
421 else if (op == Operation::add)
423 interpolation(BoundaryWrapper(model_, entity, wLocal.order(), dirichletBlocks_[global][l]), valuesModel);
424 values[ localDof ] += valuesModel[ localDof ];
427 assert( (
unsigned int)localDof < wLocal.size() );
428 wLocal[ localDof ] = values[ localDof ];
436 void updateDirichletDofs()
const
438 if( sequence_ != space_.sequence() )
441 if( ! model_.hasDirichletBoundary() )
443 hasDirichletDofs_ = false ;
448 const int blocks = space_.blockMapper().size() ;
449 dirichletBlocks_.resize( blocks );
450 for(
int i=0; i<blocks; ++i )
451 dirichletBlocks_[ i ].fill(0);
453 typedef typename DiscreteFunctionSpaceType :: IteratorType IteratorType;
454 typedef typename IteratorType :: Entity EntityType;
456 bool hasDirichletBoundary =
false;
457 const IteratorType end = space_.end();
458 for( IteratorType it = space_.begin(); it != end; ++it )
460 const EntityType &entity = *it;
462 if( entity.hasBoundaryIntersections() )
464 hasDirichletBoundary |= searchEntityDirichletDofs( entity, model_ );
469 sequence_ = space_.sequence();
470 if( space_.gridPart().comm().size() > 1 )
474 DirichletBuilder handle( *
this, space_ , space_.blockMapper() );
475 space_.gridPart().communicate
479 catch(
const Exception &e )
481 std::cerr << e << std::endl;
482 std::cerr <<
"Exception thrown in: " << __FILE__ <<
" line:" << __LINE__ << std::endl;
485 hasDirichletDofs_ = space_.gridPart().grid().comm().max( hasDirichletBoundary );
489 hasDirichletDofs_ = hasDirichletBoundary;
495 template<
class EntityType >
496 bool searchEntityDirichletDofs(
const EntityType &entity, ModelType& model )
const
498 typedef typename DiscreteFunctionSpaceType :: GridPartType GridPartType;
500 typedef typename GridPartType :: IntersectionIteratorType
501 IntersectionIteratorType;
503 const GridPartType &gridPart = space_.gridPart();
506 bool hasDirichletBoundary =
false;
509 std::vector<size_t> globalBlockDofs(space_.blockMapper().numDofs(entity));
510 space_.blockMapper().map(entity,globalBlockDofs);
512 std::vector<char> subEntityFilter(space_.blockMapper().numDofs(entity));
514 IntersectionIteratorType it = gridPart.ibegin( entity );
515 const IntersectionIteratorType endit = gridPart.iend( entity );
516 for( ; it != endit; ++it )
518 typedef typename IntersectionIteratorType :: Intersection IntersectionType;
519 const IntersectionType& intersection = *it;
522 if( intersection.boundary() )
525 DirichletBlock block;
527 std::array<int,dimRange> dblock;
529 model.init(intersection);
530 const bool isDirichletIntersection = model.isDirichletIntersection( intersection, dblock );
531 if (isDirichletIntersection)
534 const int face = intersection.indexInInside();
535 space_.blockMapper().onSubEntity(entity,face,1,subEntityFilter);
536 for(
unsigned int i=0;i<globalBlockDofs.size();++i)
538 assert( i < subEntityFilter.size());
539 if ( subEntityFilter[i]==0)
continue;
541 for(
int k = 0; k < localBlockSize; ++k)
543 int comp = subEntityFilter[i]-1+k;
544 assert( comp < dimRange );
545 dirichletBlocks_[ globalBlockDofs[ i ] ][k] = dblock[comp];
548 hasDirichletBoundary = true ;
554 return hasDirichletBoundary;
558 const DiscreteFunctionSpaceType& space_;
559 mutable DirichletBlockVector dirichletBlocks_;
560 mutable bool hasDirichletDofs_ ;
561 mutable int sequence_ ;
563 class DirichletBuilder;
566template <
class Model,
class Space >
567class DirichletConstraints< Model,Space > :: DirichletBuilder
568 :
public CommDataHandleIF< DirichletBuilder, int >
572 typedef typename SpaceType::BlockMapperType MapperType;
574 enum { nCodim = SpaceType :: GridType :: dimension + 1 };
577 typedef int DataType;
582 typedef DirichletConstraints< Model,Space > DirichletType;
583 const DirichletType &dirichlet_;
586 const MapperType &mapper_;
588 static const int blockSize = SpaceType::localBlockSize;
591 DirichletBuilder(
const DirichletType &dirichlet,
593 const MapperType &mapper )
594 : myRank_( space.gridPart().comm().rank() ),
595 mySize_( space.gridPart().comm().
size() ),
596 dirichlet_( dirichlet ),
601 bool contains (
int dim,
int codim )
const
603 return mapper_.contains( codim );
606 bool fixedSize (
int dim,
int codim )
const
612 template<
class MessageBuffer,
class Entity >
613 inline void gather ( MessageBuffer &buffer,
614 const Entity &entity )
const
616 unsigned int localBlocks = mapper_.numEntityDofs( entity );
617 std::vector<std::size_t> globalBlockDofs(localBlocks);
618 mapper_.mapEntityDofs( entity, globalBlockDofs );
619 assert(
size(entity) == globalBlockDofs.size()*blockSize );
620 for(
unsigned int localBlock = 0 ; localBlock < globalBlockDofs.size(); ++ localBlock )
622 int global = globalBlockDofs[localBlock];
623 for (
int r=0;r<blockSize;++r)
624 if (dirichlet_.dirichletBlocks_[ global ][r] )
635 template<
class MessageBuffer,
class EntityType >
636 inline void scatter ( MessageBuffer &buffer,
637 const EntityType &entity,
640 unsigned int localBlocks = mapper_.numEntityDofs( entity );
641 std::vector<std::size_t> globalBlockDofs(localBlocks);
642 mapper_.mapEntityDofs( entity, globalBlockDofs );
643 assert( n == globalBlockDofs.size()*blockSize );
644 assert( n ==
size(entity) );
645 for(
unsigned int localBlock = 0 ; localBlock < globalBlockDofs.size(); ++ localBlock )
647 int global = globalBlockDofs[localBlock];
648 for (
int r=0;r<blockSize;++r)
652 if ( !dirichlet_.dirichletBlocks_[ global ][r] && val == 1)
653 dirichlet_.dirichletBlocks_[ global ][r] =
true;
658 template<
class Entity >
659 size_t size (
const Entity &entity )
const
661 return blockSize * mapper_.numEntityDofs( entity );
Traits::ConstDofIteratorType ConstDofIteratorType
type of the const dof iterator
Definition: discretefunction.hh:628
ConstDofIteratorType dbegin() const
Obtain the constant iterator pointing to the first dof.
Definition: discretefunction.hh:761
ConstDofIteratorType dend() const
Obtain the constant iterator pointing to the last dof.
Definition: discretefunction.hh:773
DiscreteFunctionSpaceType::RangeType RangeType
type of range vector
Definition: discretefunction.hh:614
Traits::DofIteratorType DofIteratorType
type of the dof iterator
Definition: discretefunction.hh:626
Traits::FunctionSpaceType FunctionSpaceType
type of function space
Definition: discretefunctionspace.hh:194
A vector valued function space.
Definition: functionspace.hh:60
Definition: nonblockmapper.hh:284
forward declaration
Definition: discretefunction.hh:51
Default exception for dummy implementations.
Definition: exceptions.hh:263
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr std::bool_constant<((II==value)||...)> contains(std::integer_sequence< T, II... >, std::integral_constant< T, value >)
Checks whether or not a given sequence contains a value.
Definition: integersequence.hh:137
std::size_t fixedSize
The number of data items per index if it is fixed, 0 otherwise.
Definition: variablesizecommunicator.hh:264