"Fossies" - the Fresh Open Source Software Archive

Member "dune-istl-2.7.1/dune/istl/paamg/test/anisotropic.hh" (26 Nov 2020, 5701 Bytes) of package /linux/misc/dune/dune-istl-2.7.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. See also the last Fossies "Diffs" side-by-side code changes report for "anisotropic.hh": 2.6.0_vs_2.7.0.

    1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
    2 // vi: set et ts=4 sw=2 sts=2:
    3 #ifndef ANISOTROPIC_HH
    4 #define  ANISOTROPIC_HH
    5 #include <dune/common/fmatrix.hh>
    6 #include <dune/common/unused.hh>
    7 #include <dune/common/parallel/indexset.hh>
    8 #include <dune/common/parallel/plocalindex.hh>
    9 #include <dune/common/parallel/communication.hh>
   10 #include <dune/common/scalarmatrixview.hh>
   11 #include <dune/istl/bcrsmatrix.hh>
   12 #include <dune/istl/owneroverlapcopy.hh>
   13 
   14 typedef Dune::OwnerOverlapCopyAttributeSet GridAttributes;
   15 typedef GridAttributes::AttributeSet GridFlag;
   16 typedef Dune::ParallelLocalIndex<GridFlag> LocalIndex;
   17 
   18 template<class M, class G, class L, int n>
   19 void setupPattern(int N, M& mat, Dune::ParallelIndexSet<G,L,n>& indices, int overlapStart, int overlapEnd,
   20                   int start, int end);
   21 
   22 template<class M>
   23 void fillValues(int N, M& mat, int overlapStart, int overlapEnd, int start, int end);
   24 
   25 
   26 template<class M, class G, class L, int s>
   27 void setupPattern(int N, M& mat, Dune::ParallelIndexSet<G,L,s>& indices, int overlapStart, int overlapEnd,
   28                   int start, int end)
   29 {
   30   int n = overlapEnd - overlapStart;
   31 
   32   typename M::CreateIterator iter = mat.createbegin();
   33   indices.beginResize();
   34 
   35   for(int j=0; j < N; j++)
   36     for(int i=overlapStart; i < overlapEnd; i++, ++iter) {
   37       int global = j*N+i;
   38       GridFlag flag = GridAttributes::owner;
   39       bool isPublic = false;
   40 
   41       if((i<start && i > 0) || (i>= end && i < N-1))
   42         flag=GridAttributes::copy;
   43 
   44       if(i<start+1 || i>= end-1) {
   45         isPublic = true;
   46         indices.add(global, LocalIndex(iter.index(), flag, isPublic));
   47       }
   48 
   49 
   50       iter.insert(iter.index());
   51 
   52       // i direction
   53       if(i > overlapStart )
   54         // We have a left neighbour
   55         iter.insert(iter.index()-1);
   56 
   57       if(i < overlapEnd-1)
   58         // We have a rigt neighbour
   59         iter.insert(iter.index()+1);
   60 
   61       // j direction
   62       // Overlap is a dirichlet border, discard neighbours
   63       if(flag != GridAttributes::copy) {
   64         if(j>0)
   65           // lower neighbour
   66           iter.insert(iter.index()-n);
   67         if(j<N-1)
   68           // upper neighbour
   69           iter.insert(iter.index()+n);
   70       }
   71     }
   72   indices.endResize();
   73 }
   74 
   75 template<class M, class T>
   76 void fillValues(int N, M& mat, int overlapStart, int overlapEnd, int start, int end, T eps)
   77 {
   78   DUNE_UNUSED_PARAMETER(N);
   79   typedef typename M::block_type Block;
   80   Block dval(0), bone(0), bmone(0), beps(0);
   81 
   82   auto setDiagonal = [](auto&& scalarOrMatrix, const auto& value) {
   83     auto&& matrix = Dune::Impl::asMatrix(scalarOrMatrix);
   84     for (auto rowIt = matrix.begin(); rowIt != matrix.end(); ++rowIt)
   85       (*rowIt)[rowIt.index()] = value;
   86   };
   87 
   88   setDiagonal(dval, 2.0+2.0*eps);
   89   setDiagonal(bone, 1.0);
   90   setDiagonal(bmone, -1.0);
   91   setDiagonal(beps, -eps);
   92 
   93   int n = overlapEnd-overlapStart;
   94   typedef typename M::ColIterator ColIterator;
   95   typedef typename M::RowIterator RowIterator;
   96 
   97   for (RowIterator i = mat.begin(); i != mat.end(); ++i) {
   98     // calculate coordinate
   99     int y = i.index() / n;
  100     int x = overlapStart + i.index() - y * n;
  101 
  102     ColIterator endi = (*i).end();
  103 
  104     if(x<start || x >= end) { // || x==0 || x==N-1 || y==0 || y==N-1){
  105       // overlap node is dirichlet border
  106       ColIterator j = (*i).begin();
  107 
  108       for(; j.index() < i.index(); ++j)
  109         *j=0;
  110 
  111       *j = bone;
  112 
  113       for(++j; j != endi; ++j)
  114         *j=0;
  115     }else{
  116       for(ColIterator j = (*i).begin(); j != endi; ++j)
  117         if(j.index() == i.index())
  118           *j=dval;
  119         else if(j.index()+1==i.index() || j.index()-1==i.index())
  120           *j=beps;
  121         else
  122           *j=bmone;
  123     }
  124   }
  125 }
  126 
  127 template<class V, class G, class L, int s>
  128 void setBoundary(V& lhs, V& rhs, const G& n, Dune::ParallelIndexSet<G,L,s>& indices)
  129 {
  130   typedef typename Dune::ParallelIndexSet<G,L,s>::const_iterator Iter;
  131   for(Iter i=indices.begin(); i != indices.end(); ++i) {
  132     G x = i->global()/n;
  133     G y = i->global()%n;
  134 
  135     if(x==0 || y ==0 || x==n-1 || y==n-1) {
  136       //double h = 1.0 / ((double) (n-1));
  137       lhs[i->local()]=rhs[i->local()]=0; //((double)x)*((double)y)*h*h;
  138     }
  139   }
  140 }
  141 
  142 template<class V, class G>
  143 void setBoundary(V& lhs, V& rhs, const G& N)
  144 {
  145   for(int j=0; j < N; ++j)
  146     for(int i=0; i < N; i++)
  147       if(i==0 || j ==0 || i==N-1 || j==N-1)
  148         lhs[j*N+i]=rhs[j*N+i]=0;
  149 }
  150 
  151 /**
  152  * \tparam M A matrix type
  153  */
  154 template<class MatrixEntry, class G, class L, class C, int s>
  155 Dune::BCRSMatrix<MatrixEntry> setupAnisotropic2d(int N, Dune::ParallelIndexSet<G,L,s>& indices, const Dune::CollectiveCommunication<C>& p, int *nout, typename Dune::BCRSMatrix<MatrixEntry>::field_type eps=1.0)
  156 {
  157   int procs=p.size(), rank=p.rank();
  158 
  159   using BCRSMat = Dune::BCRSMatrix<MatrixEntry>;
  160 
  161   // calculate size of local matrix in the distributed direction
  162   int start, end, overlapStart, overlapEnd;
  163 
  164   int n = N/procs; // number of unknowns per process
  165   int bigger = N%procs; // number of process with n+1 unknows
  166 
  167   // Compute owner region
  168   if(rank<bigger) {
  169     start = rank*(n+1);
  170     end   = start+(n+1);
  171   }else{
  172     start = bigger*(n+1) + (rank-bigger) * n;
  173     end   = start+n;
  174   }
  175 
  176   // Compute overlap region
  177   if(start>0)
  178     overlapStart = start - 1;
  179   else
  180     overlapStart = start;
  181 
  182   if(end<N)
  183     overlapEnd = end + 1;
  184   else
  185     overlapEnd = end;
  186 
  187   int noKnown = overlapEnd-overlapStart;
  188 
  189   *nout = noKnown;
  190 
  191   BCRSMat mat(noKnown*N, noKnown*N, noKnown*N*5, BCRSMat::row_wise);
  192 
  193   setupPattern(N, mat, indices, overlapStart, overlapEnd, start, end);
  194   fillValues(N, mat, overlapStart, overlapEnd, start, end, eps);
  195 
  196   //  Dune::printmatrix(std::cout,mat,"aniso 2d","row",9,1);
  197 
  198   return mat;
  199 }
  200 #endif