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<
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);
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(
180 typename CU=EmptyTransformation,
181 typename CV=EmptyTransformation>
217 typename std::enable_if<
225 typedef typename CV::const_iterator global_col_iterator;
227 typedef typename global_col_iterator::value_type::first_type GlobalIndex;
228 const GlobalIndex & contributor = cit->first;
230 typedef typename global_col_iterator::value_type::second_type ContributedMap;
231 typedef typename ContributedMap::const_iterator global_row_iterator;
232 const ContributedMap & contributed = cit->second;
233 global_row_iterator it = contributed.begin();
234 global_row_iterator eit = contributed.end();
241 x[it->first] += it->second * x[contributor];
253 typename std::enable_if<
268 typename std::enable_if<
276 typedef typename CV::const_iterator global_col_iterator;
278 typedef typename global_col_iterator::value_type::first_type GlobalIndex;
279 const GlobalIndex & contributor = cit->first;
281 typedef typename global_col_iterator::value_type::second_type ContributedMap;
282 typedef typename ContributedMap::const_iterator global_row_iterator;
283 const ContributedMap & contributed = cit->second;
284 global_row_iterator it = contributed.begin();
285 global_row_iterator eit = contributed.end();
295 x[contributor] += it->second * x[it->first];
302 typename std::enable_if<
315 template<
typename GCView,
typename T>
318 for (
int i = 0; i < localcontainer.N(); ++i)
319 for (
int j = 0; j < localcontainer.M(); ++j)
320 localcontainer(i,j) = globalcontainer_view(i,j);
324 template<
typename T,
typename GCView>
327 for (
int i = 0; i < localcontainer.N(); ++i)
328 for (
int j = 0; j < localcontainer.M(); ++j)
329 globalcontainer_view(i,j) = localcontainer(i,j);
333 template<
typename T,
typename GCView>
336 for (
int i = 0; i < localcontainer.N(); ++i)
337 for (
int j = 0; j < localcontainer.M(); ++j)
338 globalcontainer_view.add(i,j,localcontainer(i,j));
342 template<
typename M,
typename GCView>
343 typename std::enable_if<
349 scatter_jacobian(M& local_container, GCView& global_container_view,
bool symmetric_mode)
const
351 typedef typename GCView::RowIndexCache LFSVIndexCache;
352 typedef typename GCView::ColIndexCache LFSUIndexCache;
354 const LFSVIndexCache& lfsv_indices = global_container_view.rowIndexCache();
355 const LFSUIndexCache& lfsu_indices = global_container_view.colIndexCache();
357 if (lfsv_indices.constraintsCachingEnabled() && lfsu_indices.constraintsCachingEnabled())
361 etadd(local_container,global_container_view);
365 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
366 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
368 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
369 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
374 typedef typename LFSUIndexCache::ContainerIndex CI;
376 for (
size_t j = 0; j < lfsu_indices.size(); ++j)
378 const CI& container_index = lfsu_indices.containerIndex(j);
379 const typename CU::const_iterator cit =
pconstraintsu->find(container_index);
383 assert(cit->second.empty());
385 for (
size_t i = 0; i < lfsv_indices.size(); ++i)
389 local_container(lfsv,i,lfsu,j) = 0.0;
397 for (
auto it = local_container.begin(); it != local_container.end(); ++it)
402 global_container_view.add(it.row(),it.col(),*it);
408 template<
typename M,
typename GCView>
409 typename std::enable_if<
415 scatter_jacobian(M& local_container, GCView& global_container_view,
bool symmetric_mode)
const
419 for (
auto it = local_container.begin(); it != local_container.end(); ++it)
424 global_container_view.add(it.row(),it.col(),*it);
431 template<
typename M,
typename GCView>
435 typedef typename GCView::RowIndexCache LFSVIndexCache;
436 typedef typename GCView::ColIndexCache LFSUIndexCache;
438 const LFSVIndexCache& lfsv_indices = globalcontainer_view.rowIndexCache();
439 const LFSUIndexCache& lfsu_indices = globalcontainer_view.colIndexCache();
441 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
442 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
444 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
445 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
447 for (
size_t j = 0; j < lfsu_indices.size(); ++j)
449 if (lfsu_indices.isConstrained(j) && lfsu_indices.isDirichletConstraint(j))
452 for (
size_t i = 0; i < lfsv_indices.size(); ++i)
456 localcontainer(lfsv,i,lfsu,j) = 0.0;
462 etadd(localcontainer,globalcontainer_view);
466 template<
typename M,
typename GCView>
467 void etadd (
const M& localcontainer, GCView& globalcontainer_view)
const
470 typedef typename M::value_type T;
472 typedef typename GCView::RowIndexCache LFSVIndexCache;
473 typedef typename GCView::ColIndexCache LFSUIndexCache;
475 const LFSVIndexCache& lfsv_indices = globalcontainer_view.rowIndexCache();
476 const LFSUIndexCache& lfsu_indices = globalcontainer_view.colIndexCache();
478 typedef typename LFSVIndexCache::LocalFunctionSpace LFSV;
479 const LFSV& lfsv = lfsv_indices.localFunctionSpace();
481 typedef typename LFSUIndexCache::LocalFunctionSpace LFSU;
482 const LFSU& lfsu = lfsu_indices.localFunctionSpace();
484 for (
size_t i = 0; i<lfsv_indices.size(); ++i)
485 for (
size_t j = 0; j<lfsu_indices.size(); ++j)
488 if (localcontainer(lfsv,i,lfsu,j) == 0.0)
491 const bool constrained_v = lfsv_indices.isConstrained(i);
492 const bool constrained_u = lfsu_indices.isConstrained(j);
494 typedef typename LFSVIndexCache::ConstraintsIterator VConstraintsIterator;
495 typedef typename LFSUIndexCache::ConstraintsIterator UConstraintsIterator;
499 if (lfsv_indices.isDirichletConstraint(i))
502 for (VConstraintsIterator vcit = lfsv_indices.constraintsBegin(i); vcit != lfsv_indices.constraintsEnd(i); ++vcit)
505 if (lfsu_indices.isDirichletConstraint(j))
507 T
value = localcontainer(lfsv,i,lfsu,j) * vcit->weight();
509 globalcontainer_view.add(vcit->containerIndex(),j,
value);
513 for (UConstraintsIterator ucit = lfsu_indices.constraintsBegin(j); ucit != lfsu_indices.constraintsEnd(j); ++ucit)
515 T
value = localcontainer(lfsv,i,lfsu,j) * vcit->weight() * ucit->weight();
517 globalcontainer_view.add(vcit->containerIndex(),ucit->containerIndex(),
value);
523 T
value = localcontainer(lfsv,i,lfsu,j) * vcit->weight();
525 globalcontainer_view.add(vcit->containerIndex(),j,
value);
532 if (lfsu_indices.isDirichletConstraint(j))
534 T
value = localcontainer(lfsv,i,lfsu,j);
536 globalcontainer_view.add(i,j,
value);
540 for (UConstraintsIterator ucit = lfsu_indices.constraintsBegin(j); ucit != lfsu_indices.constraintsEnd(j); ++ucit)
542 T
value = localcontainer(lfsv,i,lfsu,j) * ucit->weight();
544 globalcontainer_view.add(i,ucit->containerIndex(),
value);
549 globalcontainer_view.add(i,j,localcontainer(lfsv,i,lfsu,j));
555 template<
typename Pattern,
typename RI,
typename CI>
556 typename std::enable_if<
562 pattern.add_link(ri,ci);
565 template<
typename Pattern,
typename RI,
typename CI>
566 typename std::enable_if<
576 template<
typename P,
typename LFSVIndices,
typename LFSUIndices,
typename Index>
578 const LFSVIndices& lfsv_indices, Index i,
579 const LFSUIndices& lfsu_indices, Index j)
const
581 typedef typename LFSVIndices::ConstraintsIterator VConstraintsIterator;
582 typedef typename LFSUIndices::ConstraintsIterator UConstraintsIterator;
584 const bool constrained_v = lfsv_indices.isConstrained(i);
585 const bool constrained_u = lfsu_indices.isConstrained(j);
588 lfsv_indices.containerIndex(i),
589 lfsu_indices.containerIndex(j)
594 if (!constrained_u || lfsu_indices.isDirichletConstraint(j))
596 globalpattern.add_link(lfsv_indices.containerIndex(i),lfsu_indices.containerIndex(j));
600 for (UConstraintsIterator gurit = lfsu_indices.constraintsBegin(j); gurit != lfsu_indices.constraintsEnd(j); ++gurit)
601 globalpattern.add_link(lfsv_indices.containerIndex(i),gurit->containerIndex());
606 if (lfsv_indices.isDirichletConstraint(i))
608 globalpattern.add_link(lfsv_indices.containerIndex(i),lfsu_indices.containerIndex(j));
612 for(VConstraintsIterator gvrit = lfsv_indices.constraintsBegin(i); gvrit != lfsv_indices.constraintsEnd(i); ++gvrit)
614 if (!constrained_u || lfsu_indices.isDirichletConstraint(j))
616 globalpattern.add_link(gvrit->containerIndex(),lfsu_indices.containerIndex(j));
620 for (UConstraintsIterator gurit = lfsu_indices.constraintsBegin(j); gurit != lfsu_indices.constraintsEnd(j); ++gurit)
621 globalpattern.add_link(gvrit->containerIndex(),gurit->containerIndex());
631 template<
typename GFSV,
typename GC,
typename C>
634 typedef typename C::const_iterator global_row_iterator;
635 for (global_row_iterator cit = c.begin(); cit != c.end(); ++cit)
636 globalcontainer.clear_row(cit->first,1);
639 template<
typename GFSV,
typename GC>
644 template<
typename GFSV,
typename GC>
647 globalcontainer.flush();
649 globalcontainer.finalize();
659 template<
typename B,
typename CU,
typename CV>
661 template<
typename B,
typename CU,
typename CV>
667 #endif // DUNE_PDELAB_GRIDOPERATOR_COMMON_ASSEMBLERUTILITIES_HH