4 #ifndef DUNE_PDELAB_GRIDOPERATOR_COMMON_ASSEMBLERUTILITIES_HH
5 #define DUNE_PDELAB_GRIDOPERATOR_COMMON_ASSEMBLERUTILITIES_HH
47 typedef typename GO::Traits::Domain
Domain;
57 typedef typename GO::Traits::Range
Range;
67 typedef typename GO::Traits::Jacobian
Jacobian;
70 typedef typename MatrixBackend::template Pattern<
106 : std::tuple<int,int>(
i,
j)
110 inline int i ()
const
112 return std::get<0>(*
this);
116 inline int j ()
const
118 return std::get<1>(*
this);
124 std::get<0>(*
this) =
i;
125 std::get<1>(*
this) =
j;
136 :
public std::vector<SparsityLink>
140 using std::vector<SparsityLink>::push_back;
154 template<
typename LFSV,
typename LFSU>
155 void addLink(
const LFSV& lfsv, std::size_t i,
const LFSU& lfsu, std::size_t j)
157 std::vector<SparsityLink>::push_back(
186 typename CU=EmptyTransformation,
187 typename CV=EmptyTransformation>
223 typename std::enable_if<
231 typedef typename CV::const_iterator global_col_iterator;
233 typedef typename global_col_iterator::value_type::first_type GlobalIndex;
234 const GlobalIndex & contributor = cit->first;
236 typedef typename global_col_iterator::value_type::second_type ContributedMap;
237 typedef typename ContributedMap::const_iterator global_row_iterator;
238 const ContributedMap & contributed = cit->second;
239 global_row_iterator it = contributed.begin();
240 global_row_iterator eit = contributed.end();
247 x[it->first] += it->second * x[contributor];
259 typename std::enable_if<
274 typename std::enable_if<
282 typedef typename CV::const_iterator global_col_iterator;
284 typedef typename global_col_iterator::value_type::first_type GlobalIndex;
285 const GlobalIndex & contributor = cit->first;
287 typedef typename global_col_iterator::value_type::second_type ContributedMap;
288 typedef typename ContributedMap::const_iterator global_row_iterator;
289 const ContributedMap & contributed = cit->second;
290 global_row_iterator it = contributed.begin();
291 global_row_iterator eit = contributed.end();
301 x[contributor] += it->second * x[it->first];
308 typename std::enable_if<
321 template<
typename GCView,
typename T>
324 for (
int i = 0; i < localcontainer.N(); ++i)
325 for (
int j = 0; j < localcontainer.M(); ++j)
326 localcontainer(i,j) = globalcontainer_view(i,j);
330 template<
typename T,
typename GCView>
333 for (
int i = 0; i < localcontainer.N(); ++i)
334 for (
int j = 0; j < localcontainer.M(); ++j)
335 globalcontainer_view(i,j) = localcontainer(i,j);
339 template<
typename T,
typename GCView>
342 for (
int i = 0; i < localcontainer.N(); ++i)
343 for (
int j = 0; j < localcontainer.M(); ++j)
344 globalcontainer_view.add(i,j,localcontainer(i,j));
348 template<
typename M,
typename GCView>
349 typename std::enable_if<
355 scatter_jacobian(M& local_container, GCView& global_container_view,
bool symmetric_mode)
const
357 typedef typename GCView::RowIndexCache LFSVIndexCache;
358 typedef typename GCView::ColIndexCache LFSUIndexCache;
360 const LFSVIndexCache& lfsv_indices = global_container_view.rowIndexCache();
361 const LFSUIndexCache& lfsu_indices = global_container_view.colIndexCache();
363 if (lfsv_indices.constraintsCachingEnabled() && lfsu_indices.constraintsCachingEnabled())
367 etadd(local_container,global_container_view);
371 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
372 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
374 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
375 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
380 typedef typename LFSUIndexCache::ContainerIndex CI;
382 for (
size_t j = 0; j < lfsu_indices.size(); ++j)
384 const CI& container_index = lfsu_indices.containerIndex(j);
385 const typename CU::const_iterator cit =
pconstraintsu->find(container_index);
389 assert(cit->second.empty());
391 for (
size_t i = 0; i < lfsv_indices.size(); ++i)
395 local_container(lfsv,i,lfsu,j) = 0.0;
403 for (
auto it = local_container.begin(); it != local_container.end(); ++it)
408 global_container_view.add(it.row(),it.col(),*it);
414 template<
typename M,
typename GCView>
415 typename std::enable_if<
421 scatter_jacobian(M& local_container, GCView& global_container_view,
bool symmetric_mode)
const
425 for (
auto it = local_container.begin(); it != local_container.end(); ++it)
430 global_container_view.add(it.row(),it.col(),*it);
437 template<
typename M,
typename GCView>
441 typedef typename GCView::RowIndexCache LFSVIndexCache;
442 typedef typename GCView::ColIndexCache LFSUIndexCache;
444 const LFSVIndexCache& lfsv_indices = globalcontainer_view.rowIndexCache();
445 const LFSUIndexCache& lfsu_indices = globalcontainer_view.colIndexCache();
447 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
448 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
450 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
451 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
453 for (
size_t j = 0; j < lfsu_indices.size(); ++j)
455 if (lfsu_indices.isConstrained(j) && lfsu_indices.isDirichletConstraint(j))
458 for (
size_t i = 0; i < lfsv_indices.size(); ++i)
462 localcontainer(lfsv,i,lfsu,j) = 0.0;
468 etadd(localcontainer,globalcontainer_view);
472 template<
typename M,
typename GCView>
473 void etadd (
const M& localcontainer, GCView& globalcontainer_view)
const
476 typedef typename M::value_type T;
478 typedef typename GCView::RowIndexCache LFSVIndexCache;
479 typedef typename GCView::ColIndexCache LFSUIndexCache;
481 const LFSVIndexCache& lfsv_indices = globalcontainer_view.rowIndexCache();
482 const LFSUIndexCache& lfsu_indices = globalcontainer_view.colIndexCache();
484 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
485 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
487 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
488 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
490 for (
size_t i = 0; i<lfsv_indices.size(); ++i)
491 for (
size_t j = 0; j<lfsu_indices.size(); ++j)
494 if (localcontainer(lfsv,i,lfsu,j) == 0.0)
497 const bool constrained_v = lfsv_indices.isConstrained(i);
498 const bool constrained_u = lfsu_indices.isConstrained(j);
500 typedef typename LFSVIndexCache::ConstraintsIterator VConstraintsIterator;
501 typedef typename LFSUIndexCache::ConstraintsIterator UConstraintsIterator;
505 if (lfsv_indices.isDirichletConstraint(i))
508 for (VConstraintsIterator vcit = lfsv_indices.constraintsBegin(i); vcit != lfsv_indices.constraintsEnd(i); ++vcit)
511 if (lfsu_indices.isDirichletConstraint(j))
513 T
value = localcontainer(lfsv,i,lfsu,j) * vcit->weight();
515 globalcontainer_view.add(vcit->containerIndex(),j,
value);
519 for (UConstraintsIterator ucit = lfsu_indices.constraintsBegin(j); ucit != lfsu_indices.constraintsEnd(j); ++ucit)
521 T
value = localcontainer(lfsv,i,lfsu,j) * vcit->weight() * ucit->weight();
523 globalcontainer_view.add(vcit->containerIndex(),ucit->containerIndex(),
value);
529 T
value = localcontainer(lfsv,i,lfsu,j) * vcit->weight();
531 globalcontainer_view.add(vcit->containerIndex(),j,
value);
538 if (lfsu_indices.isDirichletConstraint(j))
540 T
value = localcontainer(lfsv,i,lfsu,j);
542 globalcontainer_view.add(i,j,
value);
546 for (UConstraintsIterator ucit = lfsu_indices.constraintsBegin(j); ucit != lfsu_indices.constraintsEnd(j); ++ucit)
548 T
value = localcontainer(lfsv,i,lfsu,j) * ucit->weight();
550 globalcontainer_view.add(i,ucit->containerIndex(),
value);
555 globalcontainer_view.add(i,j,localcontainer(lfsv,i,lfsu,j));
561 template<
typename Pattern,
typename RI,
typename CI>
562 typename std::enable_if<
568 pattern.add_link(ri,ci);
571 template<
typename Pattern,
typename RI,
typename CI>
572 typename std::enable_if<
582 template<
typename P,
typename LFSVIndices,
typename LFSUIndices,
typename Index>
584 const LFSVIndices& lfsv_indices, Index i,
585 const LFSUIndices& lfsu_indices, Index j)
const
587 typedef typename LFSVIndices::ConstraintsIterator VConstraintsIterator;
588 typedef typename LFSUIndices::ConstraintsIterator UConstraintsIterator;
590 const bool constrained_v = lfsv_indices.isConstrained(i);
591 const bool constrained_u = lfsu_indices.isConstrained(j);
594 lfsv_indices.containerIndex(i),
595 lfsu_indices.containerIndex(j)
600 if (!constrained_u || lfsu_indices.isDirichletConstraint(j))
602 globalpattern.add_link(lfsv_indices.containerIndex(i),lfsu_indices.containerIndex(j));
606 for (UConstraintsIterator gurit = lfsu_indices.constraintsBegin(j); gurit != lfsu_indices.constraintsEnd(j); ++gurit)
607 globalpattern.add_link(lfsv_indices.containerIndex(i),gurit->containerIndex());
612 if (lfsv_indices.isDirichletConstraint(i))
614 globalpattern.add_link(lfsv_indices.containerIndex(i),lfsu_indices.containerIndex(j));
618 for(VConstraintsIterator gvrit = lfsv_indices.constraintsBegin(i); gvrit != lfsv_indices.constraintsEnd(i); ++gvrit)
620 if (!constrained_u || lfsu_indices.isDirichletConstraint(j))
622 globalpattern.add_link(gvrit->containerIndex(),lfsu_indices.containerIndex(j));
626 for (UConstraintsIterator gurit = lfsu_indices.constraintsBegin(j); gurit != lfsu_indices.constraintsEnd(j); ++gurit)
627 globalpattern.add_link(gvrit->containerIndex(),gurit->containerIndex());
637 template<
typename GFSV,
typename GC,
typename C>
640 typedef typename C::const_iterator global_row_iterator;
641 for (global_row_iterator cit = c.begin(); cit != c.end(); ++cit)
642 globalcontainer.clear_row(cit->first,1);
645 template<
typename GFSV,
typename GC>
650 template<
typename GFSV,
typename GC>
653 globalcontainer.flush();
655 globalcontainer.finalize();
665 template<
typename B,
typename CU,
typename CV>
667 template<
typename B,
typename CU,
typename CV>
673 #endif // DUNE_PDELAB_GRIDOPERATOR_COMMON_ASSEMBLERUTILITIES_HH