4#ifndef DUNE_PDELAB_GRIDOPERATOR_COMMON_ASSEMBLERUTILITIES_HH
5#define DUNE_PDELAB_GRIDOPERATOR_COMMON_ASSEMBLERUTILITIES_HH
10#include <dune/pdelab/constraints/common/constraintstransformation.hh>
11#include <dune/pdelab/gridoperator/common/localmatrix.hh>
24 struct HasDoSkipEntity
27 auto require(LO&& lo) ->
decltype(
28 Concept::requireConvertible<bool>(LO::doSkipEntity)
32 struct HasDoSkipIntersection
35 auto require(LO&& lo) ->
decltype(
36 Concept::requireConvertible<bool>(LO::doSkipIntersection)
75 typedef typename GO::Traits::Domain
Domain;
85 typedef typename GO::Traits::Range
Range;
95 typedef typename GO::Traits::Jacobian
Jacobian;
98 typedef typename MatrixBackend::template Pattern<
134 :
std::tuple<int,int>(
i,
j)
138 inline int i ()
const
140 return std::get<0>(*
this);
144 inline int j ()
const
146 return std::get<1>(*
this);
152 std::get<0>(*
this) =
i;
153 std::get<1>(*
this) =
j;
164 :
public std::vector<SparsityLink>
182 template<
typename LFSV,
typename LFSU>
183 void addLink(
const LFSV& lfsv, std::size_t i,
const LFSU& lfsu, std::size_t j)
214 typename CU=EmptyTransformation,
215 typename CV=EmptyTransformation>
220 typedef typename B::size_type SizeType;
224 : pconstraintsu(&emptyconstraintsu), pconstraintsv(&emptyconstraintsv)
229 :pconstraintsu(&cu), pconstraintsv(&cv)
235 return *pconstraintsu;
241 return *pconstraintsv;
251 typename std::enable_if<
259 typedef typename CV::const_iterator global_col_iterator;
260 for (global_col_iterator cit=pconstraintsv->begin(); cit!=pconstraintsv->end(); ++cit){
261 typedef typename global_col_iterator::value_type::first_type GlobalIndex;
262 const GlobalIndex & contributor = cit->first;
264 typedef typename global_col_iterator::value_type::second_type ContributedMap;
265 typedef typename ContributedMap::const_iterator global_row_iterator;
266 const ContributedMap & contributed = cit->second;
267 global_row_iterator it = contributed.begin();
268 global_row_iterator eit = contributed.end();
275 x[it->first] += it->second * x[contributor];
280 for (global_col_iterator cit=pconstraintsv->begin(); cit!=pconstraintsv->end(); ++cit)
287 typename std::enable_if<
302 typename std::enable_if<
310 typedef typename CV::const_iterator global_col_iterator;
311 for (global_col_iterator cit=pconstraintsv->begin(); cit!=pconstraintsv->end(); ++cit){
312 typedef typename global_col_iterator::value_type::first_type GlobalIndex;
313 const GlobalIndex & contributor = cit->first;
315 typedef typename global_col_iterator::value_type::second_type ContributedMap;
316 typedef typename ContributedMap::const_iterator global_row_iterator;
317 const ContributedMap & contributed = cit->second;
318 global_row_iterator it = contributed.begin();
319 global_row_iterator eit = contributed.end();
329 x[contributor] += it->second * x[it->first];
336 typename std::enable_if<
349 template<
typename GCView,
typename T>
352 for (
int i = 0; i < localcontainer.N(); ++i)
353 for (
int j = 0; j < localcontainer.M(); ++j)
354 localcontainer(i,j) = globalcontainer_view(i,j);
358 template<
typename T,
typename GCView>
361 for (
int i = 0; i < localcontainer.N(); ++i)
362 for (
int j = 0; j < localcontainer.M(); ++j)
363 globalcontainer_view(i,j) = localcontainer(i,j);
367 template<
typename T,
typename GCView>
370 for (
int i = 0; i < localcontainer.N(); ++i)
371 for (
int j = 0; j < localcontainer.M(); ++j)
372 globalcontainer_view.add(i,j,localcontainer(i,j));
376 template<
typename M,
typename GCView>
377 typename std::enable_if<
383 scatter_jacobian(M& local_container, GCView& global_container_view,
bool symmetric_mode)
const
385 typedef typename GCView::RowIndexCache LFSVIndexCache;
386 typedef typename GCView::ColIndexCache LFSUIndexCache;
388 const LFSVIndexCache& lfsv_indices = global_container_view.rowIndexCache();
389 const LFSUIndexCache& lfsu_indices = global_container_view.colIndexCache();
391 if (lfsv_indices.constraintsCachingEnabled() && lfsu_indices.constraintsCachingEnabled())
395 etadd(local_container,global_container_view);
399 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
400 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
402 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
403 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
408 typedef typename LFSUIndexCache::ContainerIndex CI;
410 for (
size_t j = 0; j < lfsu_indices.size(); ++j)
412 const CI& container_index = lfsu_indices.containerIndex(j);
413 const typename CU::const_iterator cit = pconstraintsu->find(container_index);
414 if (cit != pconstraintsu->end())
417 assert(cit->second.empty());
419 for (
size_t i = 0; i < lfsv_indices.size(); ++i)
423 local_container(lfsv,i,lfsu,j) = 0.0;
431 for (
auto it = local_container.begin(); it != local_container.end(); ++it)
436 global_container_view.add(it.row(),it.col(),*it);
442 template<
typename M,
typename GCView>
443 typename std::enable_if<
449 scatter_jacobian(M& local_container, GCView& global_container_view,
bool symmetric_mode)
const
453 for (
auto it = local_container.begin(); it != local_container.end(); ++it)
458 global_container_view.add(it.row(),it.col(),*it);
465 template<
typename M,
typename GCView>
469 typedef typename GCView::RowIndexCache LFSVIndexCache;
470 typedef typename GCView::ColIndexCache LFSUIndexCache;
472 const LFSVIndexCache& lfsv_indices = globalcontainer_view.rowIndexCache();
473 const LFSUIndexCache& lfsu_indices = globalcontainer_view.colIndexCache();
475 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
476 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
478 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
479 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
481 for (
size_t j = 0; j < lfsu_indices.size(); ++j)
483 if (lfsu_indices.isConstrained(j) && lfsu_indices.isDirichletConstraint(j))
486 for (
size_t i = 0; i < lfsv_indices.size(); ++i)
490 localcontainer(lfsv,i,lfsu,j) = 0.0;
496 etadd(localcontainer,globalcontainer_view);
500 template<
typename M,
typename GCView>
501 void etadd (
const M& localcontainer, GCView& globalcontainer_view)
const
504 typedef typename M::value_type T;
506 typedef typename GCView::RowIndexCache LFSVIndexCache;
507 typedef typename GCView::ColIndexCache LFSUIndexCache;
509 const LFSVIndexCache& lfsv_indices = globalcontainer_view.rowIndexCache();
510 const LFSUIndexCache& lfsu_indices = globalcontainer_view.colIndexCache();
512 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
513 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
515 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
516 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
518 for (
size_t i = 0; i<lfsv_indices.size(); ++i)
519 for (
size_t j = 0; j<lfsu_indices.size(); ++j)
522 if (localcontainer(lfsv,i,lfsu,j) == 0.0)
525 const bool constrained_v = lfsv_indices.isConstrained(i);
526 const bool constrained_u = lfsu_indices.isConstrained(j);
528 typedef typename LFSVIndexCache::ConstraintsIterator VConstraintsIterator;
529 typedef typename LFSUIndexCache::ConstraintsIterator UConstraintsIterator;
533 if (lfsv_indices.isDirichletConstraint(i))
536 for (VConstraintsIterator vcit = lfsv_indices.constraintsBegin(i); vcit != lfsv_indices.constraintsEnd(i); ++vcit)
539 if (lfsu_indices.isDirichletConstraint(j))
541 T value = localcontainer(lfsv,i,lfsu,j) * vcit->weight();
543 globalcontainer_view.add(vcit->containerIndex(),j,value);
547 for (UConstraintsIterator ucit = lfsu_indices.constraintsBegin(j); ucit != lfsu_indices.constraintsEnd(j); ++ucit)
549 T value = localcontainer(lfsv,i,lfsu,j) * vcit->weight() * ucit->weight();
551 globalcontainer_view.add(vcit->containerIndex(),ucit->containerIndex(),value);
557 T value = localcontainer(lfsv,i,lfsu,j) * vcit->weight();
559 globalcontainer_view.add(vcit->containerIndex(),j,value);
566 if (lfsu_indices.isDirichletConstraint(j))
568 T value = localcontainer(lfsv,i,lfsu,j);
570 globalcontainer_view.add(i,j,value);
574 for (UConstraintsIterator ucit = lfsu_indices.constraintsBegin(j); ucit != lfsu_indices.constraintsEnd(j); ++ucit)
576 T value = localcontainer(lfsv,i,lfsu,j) * ucit->weight();
578 globalcontainer_view.add(i,ucit->containerIndex(),value);
583 globalcontainer_view.add(i,j,localcontainer(lfsv,i,lfsu,j));
589 template<
typename Pattern,
typename RI,
typename CI>
590 typename std::enable_if<
591 std::is_same<RI,CI>::value
593 add_diagonal_entry(Pattern& pattern,
const RI& ri,
const CI& ci)
const
596 pattern.add_link(ri,ci);
599 template<
typename Pattern,
typename RI,
typename CI>
600 typename std::enable_if<
601 !std::is_same<RI,CI>::value
603 add_diagonal_entry(Pattern& pattern,
const RI& ri,
const CI& ci)
const
610 template<
typename P,
typename LFSVIndices,
typename LFSUIndices,
typename Index>
612 const LFSVIndices& lfsv_indices, Index i,
613 const LFSUIndices& lfsu_indices, Index j)
const
615 typedef typename LFSVIndices::ConstraintsIterator VConstraintsIterator;
616 typedef typename LFSUIndices::ConstraintsIterator UConstraintsIterator;
618 const bool constrained_v = lfsv_indices.isConstrained(i);
619 const bool constrained_u = lfsu_indices.isConstrained(j);
621 add_diagonal_entry(globalpattern,
622 lfsv_indices.containerIndex(i),
623 lfsu_indices.containerIndex(j)
628 if (!constrained_u || lfsu_indices.isDirichletConstraint(j))
630 globalpattern.add_link(lfsv_indices.containerIndex(i),lfsu_indices.containerIndex(j));
634 for (UConstraintsIterator gurit = lfsu_indices.constraintsBegin(j); gurit != lfsu_indices.constraintsEnd(j); ++gurit)
635 globalpattern.add_link(lfsv_indices.containerIndex(i),gurit->containerIndex());
640 if (lfsv_indices.isDirichletConstraint(i))
642 globalpattern.add_link(lfsv_indices.containerIndex(i),lfsu_indices.containerIndex(j));
646 for(VConstraintsIterator gvrit = lfsv_indices.constraintsBegin(i); gvrit != lfsv_indices.constraintsEnd(i); ++gvrit)
648 if (!constrained_u || lfsu_indices.isDirichletConstraint(j))
650 globalpattern.add_link(gvrit->containerIndex(),lfsu_indices.containerIndex(j));
654 for (UConstraintsIterator gurit = lfsu_indices.constraintsBegin(j); gurit != lfsu_indices.constraintsEnd(j); ++gurit)
655 globalpattern.add_link(gvrit->containerIndex(),gurit->containerIndex());
665 template<
typename GFSV,
typename GC,
typename C>
668 typedef typename C::const_iterator global_row_iterator;
669 for (global_row_iterator cit = c.begin(); cit != c.end(); ++cit)
670 globalcontainer.clear_row(cit->first,1);
673 template<
typename GFSV,
typename GC>
674 void set_trivial_rows(
const GFSV& gfsv, GC& globalcontainer,
const EmptyTransformation& c)
const
678 template<
typename GFSV,
typename GC>
679 void handle_dirichlet_constraints(
const GFSV& gfsv, GC& globalcontainer)
const
681 globalcontainer.flush();
683 globalcontainer.finalize();
687 const CU* pconstraintsu;
688 const CV* pconstraintsv;
689 static CU emptyconstraintsu;
690 static CV emptyconstraintsv;
693 template<
typename B,
typename CU,
typename CV>
694 CU LocalAssemblerBase<B,CU,CV>::emptyconstraintsu;
695 template<
typename B,
typename CU,
typename CV>
696 CV LocalAssemblerBase<B,CU,CV>::emptyconstraintsv;
Base class for local assembler.
Definition: assemblerutilities.hh:217
std::enable_if< AlwaysTrue< M >::value &&!std::is_same< CV, EmptyTransformation >::value >::type scatter_jacobian(M &local_container, GCView &global_container_view, bool symmetric_mode) const
Scatter local jacobian to global container.
Definition: assemblerutilities.hh:383
void set_trivial_rows(const GFSV &gfsv, GC &globalcontainer, const C &c) const
insert dirichlet constraints for row and assemble T^T_U in constrained rows
Definition: assemblerutilities.hh:666
void add_entry(P &globalpattern, const LFSVIndices &lfsv_indices, Index i, const LFSUIndices &lfsu_indices, Index j) const
Adding matrix entry to pattern with respect to the constraints contributions. This assembles the entr...
Definition: assemblerutilities.hh:611
void etadd_symmetric(M &localcontainer, GCView &globalcontainer_view) const
Add local matrix to global matrix, and apply Dirichlet constraints in a symmetric fashion....
Definition: assemblerutilities.hh:466
std::enable_if< AlwaysTrue< X >::value &&!std::is_same< CV, EmptyTransformation >::value >::type backtransform(X &x, const bool prerestrict=false) const
Transforms a vector from to . If prerestrict == true then is applied instead of the full transform...
Definition: assemblerutilities.hh:308
LocalAssemblerBase(const CU &cu, const CV &cv)
construct AssemblerSpace, with constraints
Definition: assemblerutilities.hh:228
void eread(const GCView &globalcontainer_view, LocalMatrix< T > &localcontainer) const
read local stiffness matrix for entity
Definition: assemblerutilities.hh:350
std::enable_if< AlwaysTrue< X >::value &&!std::is_same< CV, EmptyTransformation >::value >::type forwardtransform(X &x, const bool postrestrict=false) const
Transforms a vector from to . If postrestrict == true then is applied instead of the full transfor...
Definition: assemblerutilities.hh:257
const CU & trialConstraints() const
get the constraints on the trial grid function space
Definition: assemblerutilities.hh:233
void eadd(const LocalMatrix< T > &localcontainer, GCView &globalcontainer_view) const
write local stiffness matrix for entity
Definition: assemblerutilities.hh:368
void ewrite(const LocalMatrix< T > &localcontainer, GCView &globalcontainer_view) const
write local stiffness matrix for entity
Definition: assemblerutilities.hh:359
LocalAssemblerBase()
construct AssemblerSpace
Definition: assemblerutilities.hh:223
const CV & testConstraints() const
get the constraints on the test grid function space
Definition: assemblerutilities.hh:239
A dense matrix for storing data associated with the degrees of freedom of a pair of LocalFunctionSpac...
Definition: localmatrix.hh:184
Layout description for a sparse linear operator.
Definition: assemblerutilities.hh:165
void addLink(const LFSV &lfsv, std::size_t i, const LFSU &lfsu, std::size_t j)
Adds a link between DOF i of lfsv and DOF j of lfsu.
Definition: assemblerutilities.hh:183
Entry in sparsity pattern.
Definition: assemblerutilities.hh:126
SparsityLink()
Standard constructor for uninitialized local index.
Definition: assemblerutilities.hh:129
void set(int i, int j)
Set both components.
Definition: assemblerutilities.hh:150
int i() const
Return first component.
Definition: assemblerutilities.hh:138
int j() const
Return second component.
Definition: assemblerutilities.hh:144
SparsityLink(int i, int j)
Initialize all components.
Definition: assemblerutilities.hh:133
constexpr HybridTreePath< T..., std::size_t > push_back(const HybridTreePath< T... > &tp, std::size_t i)
Appends a run time index to a HybridTreePath.
Definition: treepath.hh:416
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integer_sequence< T, II..., T(IN)> push_back(std::integer_sequence< T, II... >, std::integral_constant< T, IN >={})
Append an index IN to the back of the sequence.
Definition: integersequence.hh:69
template which always yields a true value
Definition: typetraits.hh:134
Definition: assemblerutilities.hh:51
GO::Traits::Range Residual
The type of the range (residual).
Definition: assemblerutilities.hh:88
MatrixBackend::template Pattern< Jacobian, TestGridFunctionSpace, TrialGridFunctionSpace > MatrixPattern
The matrix pattern.
Definition: assemblerutilities.hh:102
GO::Traits::TrialGridFunctionSpace TrialGridFunctionSpace
The trial grid function space.
Definition: assemblerutilities.hh:54
GO::Traits::Jacobian Jacobian
The type of the jacobian.
Definition: assemblerutilities.hh:95
GO::Traits::DomainField DomainField
The field type of the domain (solution).
Definition: assemblerutilities.hh:72
GO::Traits::JacobianField JacobianField
The field type of the jacobian.
Definition: assemblerutilities.hh:92
GO::Traits::TestGridFunctionSpace TestGridFunctionSpace
The test grid function space.
Definition: assemblerutilities.hh:57
GO::BorderDOFExchanger BorderDOFExchanger
The helper class to exchange data on the processor boundary.
Definition: assemblerutilities.hh:105
GO::Traits::MatrixBackend MatrixBackend
The matrix backend of the grid operator.
Definition: assemblerutilities.hh:68
GO::Traits::Domain Solution
The type of the domain (solution).
Definition: assemblerutilities.hh:78
GO::Traits::RangeField RangeField
The field type of the range (residual).
Definition: assemblerutilities.hh:82
GO::Traits::Range Range
The type of the range (residual).
Definition: assemblerutilities.hh:85
GO::Traits::Domain Domain
The type of the domain (solution).
Definition: assemblerutilities.hh:75
GO::Traits::TrialGridFunctionSpaceConstraints TrialGridFunctionSpaceConstraints
The type of the trial grid function space constraints.
Definition: assemblerutilities.hh:61
GO::Traits::TestGridFunctionSpaceConstraints TestGridFunctionSpaceConstraints
The type of the test grid function space constraints.
Definition: assemblerutilities.hh:64