4 #ifndef DUNE_PDELAB_GRIDOPERATOR_COMMON_ASSEMBLERUTILITIES_HH 5 #define DUNE_PDELAB_GRIDOPERATOR_COMMON_ASSEMBLERUTILITIES_HH 61 typedef typename GO::Traits::Jacobian
Jacobian;
64 typedef typename MatrixBackend::template Pattern<
67 TrialGridFunctionSpace
100 :
std::tuple<int,int>(i,j)
104 inline int i ()
const 106 return std::get<0>(*this);
110 inline int j ()
const 112 return std::get<1>(*this);
116 void set (
int i,
int j)
118 std::get<0>(*this) = i;
119 std::get<1>(*this) = j;
130 :
public std::vector<SparsityLink>
134 using std::vector<SparsityLink>::push_back;
148 template<
typename LFSV,
typename LFSU>
149 void addLink(
const LFSV& lfsv, std::size_t i,
const LFSU& lfsu, std::size_t j)
151 std::vector<SparsityLink>::push_back(
189 : pconstraintsu(&emptyconstraintsu), pconstraintsv(&emptyconstraintsv)
194 :pconstraintsu(&cu), pconstraintsv(&cv)
200 return *pconstraintsu;
206 return *pconstraintsv;
216 typename std::enable_if<
224 typedef typename CV::const_iterator global_col_iterator;
225 for (global_col_iterator cit=pconstraintsv->begin(); cit!=pconstraintsv->end(); ++cit){
226 typedef typename global_col_iterator::value_type::first_type GlobalIndex;
227 const GlobalIndex & contributor = cit->first;
229 typedef typename global_col_iterator::value_type::second_type ContributedMap;
230 typedef typename ContributedMap::const_iterator global_row_iterator;
231 const ContributedMap & contributed = cit->second;
232 global_row_iterator it = contributed.begin();
233 global_row_iterator eit = contributed.end();
240 x[it->first] += it->second * x[contributor];
245 for (global_col_iterator cit=pconstraintsv->begin(); cit!=pconstraintsv->end(); ++cit)
252 typename std::enable_if<
267 typename std::enable_if<
275 typedef typename CV::const_iterator global_col_iterator;
276 for (global_col_iterator cit=pconstraintsv->begin(); cit!=pconstraintsv->end(); ++cit){
277 typedef typename global_col_iterator::value_type::first_type GlobalIndex;
278 const GlobalIndex & contributor = cit->first;
280 typedef typename global_col_iterator::value_type::second_type ContributedMap;
281 typedef typename ContributedMap::const_iterator global_row_iterator;
282 const ContributedMap & contributed = cit->second;
283 global_row_iterator it = contributed.begin();
284 global_row_iterator eit = contributed.end();
294 x[contributor] += it->second * x[it->first];
301 typename std::enable_if<
314 template<
typename GCView,
typename T>
317 for (
int i = 0; i < localcontainer.N(); ++i)
318 for (
int j = 0; j < localcontainer.M(); ++j)
319 localcontainer(i,j) = globalcontainer_view(i,j);
323 template<
typename T,
typename GCView>
326 for (
int i = 0; i < localcontainer.N(); ++i)
327 for (
int j = 0; j < localcontainer.M(); ++j)
328 globalcontainer_view(i,j) = localcontainer(i,j);
332 template<
typename T,
typename GCView>
335 for (
int i = 0; i < localcontainer.N(); ++i)
336 for (
int j = 0; j < localcontainer.M(); ++j)
337 globalcontainer_view.add(i,j,localcontainer(i,j));
341 template<
typename M,
typename GCView>
342 typename std::enable_if<
348 scatter_jacobian(M& local_container, GCView& global_container_view,
bool symmetric_mode)
const 350 typedef typename GCView::RowIndexCache LFSVIndexCache;
351 typedef typename GCView::ColIndexCache LFSUIndexCache;
353 const LFSVIndexCache& lfsv_indices = global_container_view.rowIndexCache();
354 const LFSUIndexCache& lfsu_indices = global_container_view.colIndexCache();
356 if (lfsv_indices.constraintsCachingEnabled() && lfsu_indices.constraintsCachingEnabled())
358 etadd_symmetric(local_container,global_container_view);
360 etadd(local_container,global_container_view);
364 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
365 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
367 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
368 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
373 typedef typename LFSUIndexCache::ContainerIndex CI;
375 for (
size_t j = 0; j < lfsu_indices.size(); ++j)
377 const CI& container_index = lfsu_indices.containerIndex(j);
378 const typename CU::const_iterator cit = pconstraintsu->find(container_index);
379 if (cit != pconstraintsu->end())
382 assert(cit->second.empty());
384 for (
size_t i = 0; i < lfsv_indices.size(); ++i)
388 local_container(lfsv,i,lfsu,j) = 0.0;
396 for (
auto it = local_container.begin(); it != local_container.end(); ++it)
401 global_container_view.add(it.row(),it.col(),*it);
407 template<
typename M,
typename GCView>
408 typename std::enable_if<
414 scatter_jacobian(M& local_container, GCView& global_container_view,
bool symmetric_mode)
const 418 for (
auto it = local_container.begin(); it != local_container.end(); ++it)
423 global_container_view.add(it.row(),it.col(),*it);
430 template<
typename M,
typename GCView>
434 typedef typename GCView::RowIndexCache LFSVIndexCache;
435 typedef typename GCView::ColIndexCache LFSUIndexCache;
437 const LFSVIndexCache& lfsv_indices = globalcontainer_view.rowIndexCache();
438 const LFSUIndexCache& lfsu_indices = globalcontainer_view.colIndexCache();
440 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
441 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
443 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
444 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
446 for (
size_t j = 0; j < lfsu_indices.size(); ++j)
448 if (lfsu_indices.isConstrained(j) && lfsu_indices.isDirichletConstraint(j))
451 for (
size_t i = 0; i < lfsv_indices.size(); ++i)
455 localcontainer(lfsv,i,lfsu,j) = 0.0;
461 etadd(localcontainer,globalcontainer_view);
465 template<
typename M,
typename GCView>
466 void etadd (
const M& localcontainer, GCView& globalcontainer_view)
const 469 typedef typename M::value_type T;
471 typedef typename GCView::RowIndexCache LFSVIndexCache;
472 typedef typename GCView::ColIndexCache LFSUIndexCache;
474 const LFSVIndexCache& lfsv_indices = globalcontainer_view.rowIndexCache();
475 const LFSUIndexCache& lfsu_indices = globalcontainer_view.colIndexCache();
477 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
478 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
480 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
481 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
483 for (
size_t i = 0; i<lfsv_indices.size(); ++i)
484 for (
size_t j = 0; j<lfsu_indices.size(); ++j)
487 if (localcontainer(lfsv,i,lfsu,j) == 0.0)
490 const bool constrained_v = lfsv_indices.isConstrained(i);
491 const bool constrained_u = lfsu_indices.isConstrained(j);
493 typedef typename LFSVIndexCache::ConstraintsIterator VConstraintsIterator;
494 typedef typename LFSUIndexCache::ConstraintsIterator UConstraintsIterator;
498 if (lfsv_indices.isDirichletConstraint(i))
501 for (VConstraintsIterator vcit = lfsv_indices.constraintsBegin(i); vcit != lfsv_indices.constraintsEnd(i); ++vcit)
504 if (lfsu_indices.isDirichletConstraint(j))
506 T
value = localcontainer(lfsv,i,lfsu,j) * vcit->weight();
508 globalcontainer_view.add(vcit->containerIndex(),j,
value);
512 for (UConstraintsIterator ucit = lfsu_indices.constraintsBegin(j); ucit != lfsu_indices.constraintsEnd(j); ++ucit)
514 T
value = localcontainer(lfsv,i,lfsu,j) * vcit->weight() * ucit->weight();
516 globalcontainer_view.add(vcit->containerIndex(),ucit->containerIndex(),
value);
522 T
value = localcontainer(lfsv,i,lfsu,j) * vcit->weight();
524 globalcontainer_view.add(vcit->containerIndex(),j,
value);
531 if (lfsu_indices.isDirichletConstraint(j))
533 T
value = localcontainer(lfsv,i,lfsu,j);
535 globalcontainer_view.add(i,j,value);
539 for (UConstraintsIterator ucit = lfsu_indices.constraintsBegin(j); ucit != lfsu_indices.constraintsEnd(j); ++ucit)
541 T
value = localcontainer(lfsv,i,lfsu,j) * ucit->weight();
543 globalcontainer_view.add(i,ucit->containerIndex(),
value);
548 globalcontainer_view.add(i,j,localcontainer(lfsv,i,lfsu,j));
554 template<
typename Pattern,
typename RI,
typename CI>
555 typename std::enable_if<
561 pattern.add_link(ri,ci);
564 template<
typename Pattern,
typename RI,
typename CI>
565 typename std::enable_if<
575 template<
typename P,
typename LFSVIndices,
typename LFSUIndices,
typename Index>
577 const LFSVIndices& lfsv_indices, Index i,
578 const LFSUIndices& lfsu_indices, Index j)
const 580 typedef typename LFSVIndices::ConstraintsIterator VConstraintsIterator;
581 typedef typename LFSUIndices::ConstraintsIterator UConstraintsIterator;
583 const bool constrained_v = lfsv_indices.isConstrained(i);
584 const bool constrained_u = lfsu_indices.isConstrained(j);
586 add_diagonal_entry(globalpattern,
587 lfsv_indices.containerIndex(i),
588 lfsu_indices.containerIndex(j)
593 if (!constrained_u || lfsu_indices.isDirichletConstraint(j))
595 globalpattern.add_link(lfsv_indices.containerIndex(i),lfsu_indices.containerIndex(j));
599 for (UConstraintsIterator gurit = lfsu_indices.constraintsBegin(j); gurit != lfsu_indices.constraintsEnd(j); ++gurit)
600 globalpattern.add_link(lfsv_indices.containerIndex(i),gurit->containerIndex());
605 if (lfsv_indices.isDirichletConstraint(i))
607 globalpattern.add_link(lfsv_indices.containerIndex(i),lfsu_indices.containerIndex(j));
611 for(VConstraintsIterator gvrit = lfsv_indices.constraintsBegin(i); gvrit != lfsv_indices.constraintsEnd(i); ++gvrit)
613 if (!constrained_u || lfsu_indices.isDirichletConstraint(j))
615 globalpattern.add_link(gvrit->containerIndex(),lfsu_indices.containerIndex(j));
619 for (UConstraintsIterator gurit = lfsu_indices.constraintsBegin(j); gurit != lfsu_indices.constraintsEnd(j); ++gurit)
620 globalpattern.add_link(gvrit->containerIndex(),gurit->containerIndex());
630 template<
typename GFSV,
typename GC,
typename C>
633 typedef typename C::const_iterator global_row_iterator;
634 for (global_row_iterator cit = c.begin(); cit != c.end(); ++cit)
635 globalcontainer.clear_row(cit->first,1);
638 template<
typename GFSV,
typename GC>
643 template<
typename GFSV,
typename GC>
646 globalcontainer.flush();
647 set_trivial_rows(gfsv,globalcontainer,*pconstraintsv);
648 globalcontainer.finalize();
658 template<
typename B,
typename CU,
typename CV>
660 template<
typename B,
typename CU,
typename CV>
666 #endif // DUNE_PDELAB_GRIDOPERATOR_COMMON_ASSEMBLERUTILITIES_HH GO::Traits::TrialGridFunctionSpace TrialGridFunctionSpace
The trial grid function space.
Definition: assemblerutilities.hh:26
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:273
std::enable_if< !std::is_same< RI, CI >::value >::type add_diagonal_entry(Pattern &pattern, const RI &ri, const CI &ci) const
Definition: assemblerutilities.hh:568
B::size_type SizeType
Definition: assemblerutilities.hh:185
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:631
const CV & testConstraints() const
get the constraints on the test grid function space
Definition: assemblerutilities.hh:204
LocalAssemblerBase(const CU &cu, const CV &cv)
construct AssemblerSpace, with constraints
Definition: assemblerutilities.hh:193
GO::Traits::MatrixBackend MatrixBackend
The matrix backend of the grid operator.
Definition: assemblerutilities.hh:40
void set_trivial_rows(const GFSV &gfsv, GC &globalcontainer, const EmptyTransformation &c) const
Definition: assemblerutilities.hh:639
GO::Traits::RangeField RangeField
The field type of the range (residual).
Definition: assemblerutilities.hh:51
void eadd(const LocalMatrix< T > &localcontainer, GCView &globalcontainer_view) const
write local stiffness matrix for entity
Definition: assemblerutilities.hh:333
GO::Traits::DomainField DomainField
The field type of the domain (solution).
Definition: assemblerutilities.hh:44
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:348
SparsityLink(int i, int j)
Initialize all components.
Definition: assemblerutilities.hh:99
Definition: assemblerutilities.hh:22
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:222
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
GO::Traits::JacobianField JacobianField
The field type of the jacobian.
Definition: assemblerutilities.hh:58
GO::Traits::TestGridFunctionSpaceConstraints TestGridFunctionSpaceConstraints
The type of the test grid function space constraints.
Definition: assemblerutilities.hh:36
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:149
Layout description for a sparse linear operator.
Definition: assemblerutilities.hh:129
int j() const
Return second component.
Definition: assemblerutilities.hh:110
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:431
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:576
GO::Traits::Jacobian Jacobian
The type of the jacobian.
Definition: assemblerutilities.hh:61
int i() const
Return first component.
Definition: assemblerutilities.hh:104
A dense matrix for storing data associated with the degrees of freedom of a pair of LocalFunctionSpac...
Definition: localmatrix.hh:184
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
const CU * pconstraintsu
Definition: assemblerutilities.hh:652
std::enable_if< std::is_same< RI, CI >::value >::type add_diagonal_entry(Pattern &pattern, const RI &ri, const CI &ci) const
Definition: assemblerutilities.hh:558
std::enable_if< AlwaysTrue< X >::value &&std::is_same< CV, EmptyTransformation >::value >::type backtransform(X &x, const bool prerestrict=false) const
Definition: assemblerutilities.hh:307
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
Definition: assemblerutilities.hh:414
GO::Traits::Range Residual
The type of the range (residual).
Definition: assemblerutilities.hh:54
GO::Traits::Domain Solution
The type of the domain (solution).
Definition: assemblerutilities.hh:47
Entry in sparsity pattern.
Definition: assemblerutilities.hh:91
void etadd(const M &localcontainer, GCView &globalcontainer_view) const
Definition: assemblerutilities.hh:466
MatrixBackend::template Pattern< Jacobian, TestGridFunctionSpace, TrialGridFunctionSpace > MatrixPattern
The matrix pattern.
Definition: assemblerutilities.hh:68
GO::BorderDOFExchanger BorderDOFExchanger
The helper class to exchange data on the processor boundary.
Definition: assemblerutilities.hh:71
void handle_dirichlet_constraints(const GFSV &gfsv, GC &globalcontainer) const
Definition: assemblerutilities.hh:644
const CU & trialConstraints() const
get the constraints on the trial grid function space
Definition: assemblerutilities.hh:198
GO::Traits::TrialGridFunctionSpaceConstraints TrialGridFunctionSpaceConstraints
The type of the trial grid function space constraints.
Definition: assemblerutilities.hh:33
static CU emptyconstraintsu
Definition: assemblerutilities.hh:654
std::enable_if< AlwaysTrue< X >::value &&std::is_same< CV, EmptyTransformation >::value >::type forwardtransform(X &x, const bool postrestrict=false) const
Definition: assemblerutilities.hh:258
Base class for local assembler.
Definition: assemblerutilities.hh:182
SparsityLink()
Standard constructor for uninitialized local index.
Definition: assemblerutilities.hh:95
const CV * pconstraintsv
Definition: assemblerutilities.hh:653
void ewrite(const LocalMatrix< T > &localcontainer, GCView &globalcontainer_view) const
write local stiffness matrix for entity
Definition: assemblerutilities.hh:324
LocalAssemblerBase()
construct AssemblerSpace
Definition: assemblerutilities.hh:188
Definition: constraintstransformation.hh:111
static CV emptyconstraintsv
Definition: assemblerutilities.hh:655
void eread(const GCView &globalcontainer_view, LocalMatrix< T > &localcontainer) const
read local stiffness matrix for entity
Definition: assemblerutilities.hh:315
GO::Traits::TestGridFunctionSpace TestGridFunctionSpace
The test grid function space.
Definition: assemblerutilities.hh:29