4 #ifndef DUNE_PDELAB_ASSEMBLERUTILITIES_HH
5 #define DUNE_PDELAB_ASSEMBLERUTILITIES_HH
61 typedef typename GO::Traits::Jacobian
Jacobian;
64 typedef typename MatrixBackend::template Pattern<
67 TrialGridFunctionSpace
93 template<
typename Intersection>
94 static Type get(
const Intersection& is)
96 return static_cast<Type>(1*is.neighbor() + 2*is.boundary());
128 :
std::tuple<int,int>(i,j)
132 inline int i ()
const
134 return std::get<0>(*this);
138 inline int j ()
const
140 return std::get<1>(*this);
146 std::get<0>(*this) =
i;
147 std::get<1>(*this) =
j;
158 :
public std::vector<SparsityLink>
162 using std::vector<SparsityLink>::push_back;
176 template<
typename LFSV,
typename LFSU>
177 void addLink(
const LFSV& lfsv, std::size_t i,
const LFSU& lfsu, std::size_t j)
179 std::vector<SparsityLink>::push_back(
208 typename CU=EmptyTransformation,
209 typename CV=EmptyTransformation>
252 typedef typename CV::const_iterator global_col_iterator;
254 typedef typename global_col_iterator::value_type::first_type GlobalIndex;
255 const GlobalIndex & contributor = cit->first;
257 typedef typename global_col_iterator::value_type::second_type ContributedMap;
258 typedef typename ContributedMap::const_iterator global_row_iterator;
259 const ContributedMap & contributed = cit->second;
260 global_row_iterator it = contributed.begin();
261 global_row_iterator eit = contributed.end();
268 x[it->first] += it->second * x[contributor];
303 typedef typename CV::const_iterator global_col_iterator;
305 typedef typename global_col_iterator::value_type::first_type GlobalIndex;
306 const GlobalIndex & contributor = cit->first;
308 typedef typename global_col_iterator::value_type::second_type ContributedMap;
309 typedef typename ContributedMap::const_iterator global_row_iterator;
310 const ContributedMap & contributed = cit->second;
311 global_row_iterator it = contributed.begin();
312 global_row_iterator eit = contributed.end();
322 x[contributor] += it->second * x[it->first];
342 template<
typename GCView,
typename T>
345 for (
int i = 0; i < localcontainer.N(); ++i)
346 for (
int j = 0; j < localcontainer.M(); ++j)
347 localcontainer(i,j) = globalcontainer_view(i,j);
351 template<
typename T,
typename GCView>
354 for (
int i = 0; i < localcontainer.N(); ++i)
355 for (
int j = 0; j < localcontainer.M(); ++j)
356 globalcontainer_view(i,j) = localcontainer(i,j);
360 template<
typename T,
typename GCView>
363 for (
int i = 0; i < localcontainer.N(); ++i)
364 for (
int j = 0; j < localcontainer.M(); ++j)
365 globalcontainer_view.add(i,j,localcontainer(i,j));
369 template<
typename M,
typename GCView>
376 scatter_jacobian(M& local_container, GCView& global_container_view,
bool symmetric_mode)
const
378 typedef typename GCView::RowIndexCache LFSVIndexCache;
379 typedef typename GCView::ColIndexCache LFSUIndexCache;
381 const LFSVIndexCache& lfsv_indices = global_container_view.rowIndexCache();
382 const LFSUIndexCache& lfsu_indices = global_container_view.colIndexCache();
384 if (lfsv_indices.constraintsCachingEnabled() && lfsu_indices.constraintsCachingEnabled())
388 etadd(local_container,global_container_view);
392 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
393 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
395 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
396 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
401 typedef typename LFSUIndexCache::ContainerIndex CI;
403 for (
size_t j = 0; j < lfsu_indices.size(); ++j)
405 const CI& container_index = lfsu_indices.containerIndex(j);
406 const typename CU::const_iterator cit =
pconstraintsu->find(container_index);
410 assert(cit->second.empty());
412 for (
size_t i = 0; i < lfsv_indices.size(); ++i)
416 local_container(lfsv,i,lfsu,j) = 0.0;
424 for (
auto it = local_container.begin(); it != local_container.end(); ++it)
429 global_container_view.add(it.row(),it.col(),*it);
435 template<
typename M,
typename GCView>
442 scatter_jacobian(M& local_container, GCView& global_container_view,
bool symmetric_mode)
const
446 for (
auto it = local_container.begin(); it != local_container.end(); ++it)
451 global_container_view.add(it.row(),it.col(),*it);
458 template<
typename M,
typename GCView>
462 typedef typename GCView::RowIndexCache LFSVIndexCache;
463 typedef typename GCView::ColIndexCache LFSUIndexCache;
465 const LFSVIndexCache& lfsv_indices = globalcontainer_view.rowIndexCache();
466 const LFSUIndexCache& lfsu_indices = globalcontainer_view.colIndexCache();
468 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
469 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
471 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
472 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
474 for (
size_t j = 0; j < lfsu_indices.size(); ++j)
476 if (lfsu_indices.isConstrained(j) && lfsu_indices.isDirichletConstraint(j))
479 for (
size_t i = 0; i < lfsv_indices.size(); ++i)
483 localcontainer(lfsv,i,lfsu,j) = 0.0;
489 etadd(localcontainer,globalcontainer_view);
493 template<
typename M,
typename GCView>
494 void etadd (
const M& localcontainer, GCView& globalcontainer_view)
const
497 typedef typename M::value_type T;
499 typedef typename GCView::RowIndexCache LFSVIndexCache;
500 typedef typename GCView::ColIndexCache LFSUIndexCache;
502 const LFSVIndexCache& lfsv_indices = globalcontainer_view.rowIndexCache();
503 const LFSUIndexCache& lfsu_indices = globalcontainer_view.colIndexCache();
505 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
506 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
508 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
509 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
511 for (
size_t i = 0; i<lfsv_indices.size(); ++i)
512 for (
size_t j = 0; j<lfsu_indices.size(); ++j)
515 if (localcontainer(lfsv,i,lfsu,j) == 0.0)
518 const bool constrained_v = lfsv_indices.isConstrained(i);
519 const bool constrained_u = lfsu_indices.isConstrained(j);
521 typedef typename LFSVIndexCache::ConstraintsIterator VConstraintsIterator;
522 typedef typename LFSUIndexCache::ConstraintsIterator UConstraintsIterator;
526 if (lfsv_indices.isDirichletConstraint(i))
529 for (VConstraintsIterator vcit = lfsv_indices.constraintsBegin(i); vcit != lfsv_indices.constraintsEnd(i); ++vcit)
532 if (lfsu_indices.isDirichletConstraint(j))
534 T
value = localcontainer(lfsv,i,lfsu,j) * vcit->weight();
536 globalcontainer_view.add(vcit->containerIndex(),j,
value);
540 for (UConstraintsIterator ucit = lfsu_indices.constraintsBegin(j); ucit != lfsu_indices.constraintsEnd(j); ++ucit)
542 T
value = localcontainer(lfsv,i,lfsu,j) * vcit->weight() * ucit->weight();
544 globalcontainer_view.add(vcit->containerIndex(),ucit->containerIndex(),
value);
550 T
value = localcontainer(lfsv,i,lfsu,j) * vcit->weight();
552 globalcontainer_view.add(vcit->containerIndex(),j,
value);
559 if (lfsu_indices.isDirichletConstraint(j))
561 T
value = localcontainer(lfsv,i,lfsu,j);
563 globalcontainer_view.add(i,j,value);
567 for (UConstraintsIterator ucit = lfsu_indices.constraintsBegin(j); ucit != lfsu_indices.constraintsEnd(j); ++ucit)
569 T
value = localcontainer(lfsv,i,lfsu,j) * ucit->weight();
571 globalcontainer_view.add(i,ucit->containerIndex(),
value);
576 globalcontainer_view.add(i,j,localcontainer(lfsv,i,lfsu,j));
582 template<
typename Pattern,
typename RI,
typename CI>
589 pattern.add_link(ri,ci);
592 template<
typename Pattern,
typename RI,
typename CI>
603 template<
typename P,
typename LFSVIndices,
typename LFSUIndices,
typename Index>
605 const LFSVIndices& lfsv_indices, Index i,
606 const LFSUIndices& lfsu_indices, Index j)
const
608 typedef typename LFSVIndices::ConstraintsIterator VConstraintsIterator;
609 typedef typename LFSUIndices::ConstraintsIterator UConstraintsIterator;
611 const bool constrained_v = lfsv_indices.isConstrained(i);
612 const bool constrained_u = lfsu_indices.isConstrained(j);
615 lfsv_indices.containerIndex(i),
616 lfsu_indices.containerIndex(j)
621 if (!constrained_u || lfsu_indices.isDirichletConstraint(j))
623 globalpattern.add_link(lfsv_indices.containerIndex(i),lfsu_indices.containerIndex(j));
627 for (UConstraintsIterator gurit = lfsu_indices.constraintsBegin(j); gurit != lfsu_indices.constraintsEnd(j); ++gurit)
628 globalpattern.add_link(lfsv_indices.containerIndex(i),gurit->containerIndex());
633 if (lfsv_indices.isDirichletConstraint(i))
635 globalpattern.add_link(lfsv_indices.containerIndex(i),lfsu_indices.containerIndex(j));
639 for(VConstraintsIterator gvrit = lfsv_indices.constraintsBegin(i); gvrit != lfsv_indices.constraintsEnd(i); ++gvrit)
641 if (!constrained_u || lfsu_indices.isDirichletConstraint(j))
643 globalpattern.add_link(gvrit->containerIndex(),lfsu_indices.containerIndex(j));
647 for (UConstraintsIterator gurit = lfsu_indices.constraintsBegin(j); gurit != lfsu_indices.constraintsEnd(j); ++gurit)
648 globalpattern.add_link(gvrit->containerIndex(),gurit->containerIndex());
658 template<
typename GFSV,
typename GC,
typename C>
661 typedef typename C::const_iterator global_row_iterator;
662 for (global_row_iterator cit = c.begin(); cit != c.end(); ++cit)
663 globalcontainer.clear_row(cit->first,1);
666 template<
typename GFSV,
typename GC>
671 template<
typename GFSV,
typename GC>
674 globalcontainer.flush();
676 globalcontainer.finalize();
686 template<
typename B,
typename CU,
typename CV>
688 template<
typename B,
typename CU,
typename CV>
694 #endif //DUNE_PDELAB_ASSEMBLERUTILITIES_HH
GO::Traits::Range Residual
The type of the range (residual).
Definition: assemblerutilities.hh:54
GO::Traits::RangeField RangeField
The field type of the range (residual).
Definition: assemblerutilities.hh:51
SparsityLink()
Standard constructor for uninitialized local index.
Definition: assemblerutilities.hh:123
const CV * pconstraintsv
Definition: assemblerutilities.hh:681
void eadd(const LocalMatrix< T > &localcontainer, GCView &globalcontainer_view) const
write local stiffness matrix for entity
Definition: assemblerutilities.hh:361
enable_if< AlwaysTrue< X >::value &&!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:250
B::size_type SizeType
Definition: assemblerutilities.hh:213
GO::Traits::Jacobian Jacobian
The type of the jacobian.
Definition: assemblerutilities.hh:61
domain boundary intersection (neighbor() == false && boundary() == true)
Definition: assemblerutilities.hh:88
GO::Traits::TrialGridFunctionSpace TrialGridFunctionSpace
The trial grid function space.
Definition: assemblerutilities.hh:26
MatrixBackend::template Pattern< Jacobian, TestGridFunctionSpace, TrialGridFunctionSpace > MatrixPattern
The matrix pattern.
Definition: assemblerutilities.hh:68
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:604
int i() const
Return first component.
Definition: assemblerutilities.hh:132
SparsityLink(int i, int j)
Initialize all components.
Definition: assemblerutilities.hh:127
enable_if< is_same< RI, CI >::value >::type add_diagonal_entry(Pattern &pattern, const RI &ri, const CI &ci) const
Definition: assemblerutilities.hh:586
periodic boundary intersection (neighbor() == true && boundary() == true)
Definition: assemblerutilities.hh:89
void etadd(const M &localcontainer, GCView &globalcontainer_view) const
Definition: assemblerutilities.hh:494
GO::Traits::DomainField DomainField
The field type of the domain (solution).
Definition: assemblerutilities.hh:44
A dense matrix for storing data associated with the degrees of freedom of a pair of LocalFunctionSpac...
Definition: localmatrix.hh:184
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:659
int j() const
Return second component.
Definition: assemblerutilities.hh:138
Entry in sparsity pattern.
Definition: assemblerutilities.hh:119
Layout description for a sparse linear operator.
Definition: assemblerutilities.hh:157
enable_if< AlwaysTrue< M >::value &&is_same< CV, EmptyTransformation >::value >::type scatter_jacobian(M &local_container, GCView &global_container_view, bool symmetric_mode) const
Definition: assemblerutilities.hh:442
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
void set(int i, int j)
Set both components.
Definition: assemblerutilities.hh:144
GO::Traits::TrialGridFunctionSpaceConstraints TrialGridFunctionSpaceConstraints
The type of the trial grid function space constraints.
Definition: assemblerutilities.hh:33
const CU * pconstraintsu
Definition: assemblerutilities.hh:680
enable_if< AlwaysTrue< X >::value &&is_same< CV, EmptyTransformation >::value >::type forwardtransform(X &x, const bool postrestrict=false) const
Definition: assemblerutilities.hh:286
void ewrite(const LocalMatrix< T > &localcontainer, GCView &globalcontainer_view) const
write local stiffness matrix for entity
Definition: assemblerutilities.hh:352
Translation helper from intersection method return values to intersection type.
Definition: assemblerutilities.hh:82
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:459
enable_if< AlwaysTrue< M >::value &&!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:376
void set_trivial_rows(const GFSV &gfsv, GC &globalcontainer, const EmptyTransformation &c) const
Definition: assemblerutilities.hh:667
void handle_dirichlet_constraints(const GFSV &gfsv, GC &globalcontainer) const
Definition: assemblerutilities.hh:672
Base class for local assembler.
Definition: assemblerutilities.hh:210
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:177
const CV & testConstraints() const
get the constraints on the test grid function space
Definition: assemblerutilities.hh:232
LocalAssemblerBase()
construct AssemblerSpace
Definition: assemblerutilities.hh:216
enable_if< AlwaysTrue< X >::value &&!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:301
Definition: adaptivity.hh:27
Definition: constraintstransformation.hh:111
void eread(const GCView &globalcontainer_view, LocalMatrix< T > &localcontainer) const
read local stiffness matrix for entity
Definition: assemblerutilities.hh:343
Definition: assemblerutilities.hh:22
LocalAssemblerBase(const CU &cu, const CV &cv)
construct AssemblerSpace, with constraints
Definition: assemblerutilities.hh:221
Type
Definition: assemblerutilities.hh:85
GO::Traits::TestGridFunctionSpaceConstraints TestGridFunctionSpaceConstraints
The type of the test grid function space constraints.
Definition: assemblerutilities.hh:36
GO::Traits::MatrixBackend MatrixBackend
The matrix backend of the grid operator.
Definition: assemblerutilities.hh:40
GO::BorderDOFExchanger BorderDOFExchanger
The helper class to exchange data on the processor boundary.
Definition: assemblerutilities.hh:71
GO::Traits::JacobianField JacobianField
The field type of the jacobian.
Definition: assemblerutilities.hh:58
enable_if< AlwaysTrue< X >::value &&is_same< CV, EmptyTransformation >::value >::type backtransform(X &x, const bool prerestrict=false) const
Definition: assemblerutilities.hh:335
skeleton intersection (neighbor() == true && boundary() == false)
Definition: assemblerutilities.hh:87
static CU emptyconstraintsu
Definition: assemblerutilities.hh:682
GO::Traits::Domain Solution
The type of the domain (solution).
Definition: assemblerutilities.hh:47
processor boundary intersection (neighbor() == false && boundary() == false)
Definition: assemblerutilities.hh:86
enable_if< !is_same< RI, CI >::value >::type add_diagonal_entry(Pattern &pattern, const RI &ri, const CI &ci) const
Definition: assemblerutilities.hh:596
const CU & trialConstraints() const
get the constraints on the trial grid function space
Definition: assemblerutilities.hh:226
GO::Traits::TestGridFunctionSpace TestGridFunctionSpace
The test grid function space.
Definition: assemblerutilities.hh:29
static CV emptyconstraintsv
Definition: assemblerutilities.hh:683