[Dune] how to use geometryInFather method
Oliver Sander
sander at mi.fu-berlin.de
Tue Aug 8 18:27:30 CEST 2006
Hi Sreejith
The expression
eit->geometryInFather().global(xi)
is correct. Now what does it mean? eit is a fine-grid element. Its
vertices have position somewhere in the 'world' of the grid. The method
geometryInFather() yields an element which is identical to eit. However,
its vertices now take positions within its father's reference element.
Lets say you have a coarse grid quad element with coorners
(1,1),(2,1),(1,2),(2,2). Subdivide it, and consider the lower left of
its four sons. Its vertices have the coordinates
(1,1),(1.5,1),(1,1.5),(1.5,1.5). Now call the method geometryInFather().
You get the embedding in the father, that is, the vertex coordinates
are now (0,0),(.5,0),(0,0.5),(0.5,0.5) (because the father's
reference element is the unit square.) This unit square constitutes
the local coordinate system for the father, but it is the world
coordinate system for the element return by geometryInFather (because
geometryInFather is constructed that way). Now in order to go
from local coordinates on the son to local coordinates on the
father (== world coordinates of geometryInFather()) you have to
call the method global().
Hope this helps,
Oliver
************************************************************************
* Oliver Sander ** email: sander at mi.fu-berlin.de *
* Freie Universität Berlin ** phone: + 49 (30) 838 75217 *
* Institut für Mathematik II ** URL : page.mi.fu-berlin.de/~sander *
* Arnimallee 6 ** -------------------------------------*
* 14195 Berlin, Germany ** Member of MATHEON (www.matheon.de) *
************************************************************************
On Tue, 8 Aug 2006, Sreejith Pulloor Kuttanikkad wrote:
> Dear Oliver,
> Thanks a lot. I think my problem is rather simpler. but i am confused a
> bit. Let me explain with a code:
> you may read the comments on the code.
>
> //---------cut here--------//
> #ifdef HAVE_CONFIG_H
> # include "config.h"
> #endif
> #include <iostream>
> #include"dune/grid/sgrid.hh" // load sgrid definition
> #include"dune/grid/io/file/vtk/vtkwriter.hh"
> int main()
> {
> // ------S Grid 2D----------
> static const int dim=2;
> typedef Dune::SGrid<dim,dim> Grid;
> int ELEMENTS[dim];double L[dim];
> //grid dimension
> L[0]=5.0; //Length
> L[1]=5.0; //width
> //discretization
> ELEMENTS[0]=1; //elements in X dir
> ELEMENTS[1]=1; //elements in Y dir
>
> //creating a father element (5 unit * 5 unit)
> Grid grid(ELEMENTS,L);
>
> //now i refine the grid uniformily
> grid.globalRefine(1);
>
> // let the point P locate at xi with respect to the element which i got
> // by subdividing the father element
>
> Dune::FieldVector<double,2> xi;
> xi[0]=0.50;
> xi[1]=0.50;
>
> int level=grid.maxLevel();
> typedef Grid::Codim<0>::LevelIterator ElementIterator;
> ElementIterator eit = grid.lbegin<0>(level);
> ElementIterator eitend = grid.lend<0>(level);
>
> // for (; eit != eitend; ++eit) // consider the very first element here
>
> // assumed P is located at xi on this element
>
>
> {
>
> std::cout<<"local coord of P with respect to father :"<<eit->geometryInFather().global(xi)<<std::endl;
>
> }
>
> // as you can see in the above, infact the above method geometryInFather().global(xi)
> // gives us the local coordinate of P
> // with respect to father geometry. but i am confused by the use of
> // global method there.
>
> }
>
> //---------cut here--------//
>
>
>
> So i simply need to get the local coordinate of the given point P with
> respect to father, when i am only given the local coordinate of P with
> respect to son.
>
>
> Thank you,
> Sreejith
>
>
> On Tue, Aug 08, 2006 at 05:47:42PM +0200, Oliver Sander wrote:
>> Hi Sreejith!
>> Das klingt als wolltest Du eine Mehrgitterrestriktion programmieren?
>> Sowas habe ich schon gemacht. Die Klasse habe ich Dir angeh?ngt.
>> Falls Du doch etwas anderes brauchst enth?lt der Code auch die
>> genaue Antwort auf Deine Frage.
>>
>> Gr??e,
>> Oliver
>>
>> ************************************************************************
>> * Oliver Sander ** email: sander at mi.fu-berlin.de *
>> * Freie Universit?t Berlin ** phone: + 49 (30) 838 75217 *
>> * Institut f?r Mathematik II ** URL : page.mi.fu-berlin.de/~sander *
>> * Arnimallee 6 ** -------------------------------------*
>> * 14195 Berlin, Germany ** Member of MATHEON (www.matheon.de) *
>> ************************************************************************
>>
>> On Tue, 8 Aug 2006, Sreejith Pulloor Kuttanikkad wrote:
>>
>>> Dear All,
>>>
>>> I would like to know how to use the geometryInFather method and in
>>> particular I want to do the following:
>>>
>>> A father element is subdivided and i am now on my element which i got from
>>> the subdivision. now i have a point "p" on this/my element say at "xi"
>>> (local coordinate with respect to this element).
>>> now i want to find the local coordinate of the point "p" with respect to
>>> the father element from "xi". how do i do that.?
>>>
>>> The documentation didnt help me much.
>>>
>>> thanks,
>>> Sreejith
>>>
>>>
>>>
>>> --
>>> Sreejith P. Kuttanikkad
>>> IWR, University of Heidelberg
>>> Room:009, Im Neuenheimer Feld-348
>>> 69120 Heidelberg,Germany.
>>> Ph :+49-6221-54-5689/4412(Office)
>>> :+49-(0)17624228904(Mob)
>>> http://hal.iwr.uni-heidelberg.de/~sreejith/
>>> -------------------------------------------
>>>
>>> _______________________________________________
>>> Dune mailing list
>>> Dune at dune-project.org
>>> http://www.dune-project.org/cgi-bin/mailman/listinfo/dune
>>>
>
>> #ifndef DUNE_MULTIGRID_TRANSFER_HH
>> #define DUNE_MULTIGRID_TRANSFER_HH
>>
>> #include <dune/istl/bcrsmatrix.hh>
>> #include <dune/common/fmatrix.hh>
>> #include <dune/common/bitfield.hh>
>>
>> namespace Dune {
>>
>>
>> /** \brief Restriction and prolongation operator for standard multigrid
>> *
>> * This class implements the standard prolongation and restriction
>> * operators for standard multigrid solvers. Restriction and prolongation
>> * of block vectors is provided. Internally, the operator is stored
>> * as a BCRSMatrix. Therefore, the template parameter DiscFuncType
>> * has to comply with the ISTL requirements.
>>
>> * \todo Currently only works for first-order Lagrangian elements!
>> */
>> template<class DiscFuncType>
>> class MultiGridTransfer {
>>
>> enum {blocksize = DiscFuncType::block_type::size};
>>
>> /** \todo Should we extract the value type from DiscFuncType? */
>> typedef FieldMatrix<double, blocksize, blocksize> MatrixBlock;
>>
>>
>> public:
>>
>> typedef BCRSMatrix<MatrixBlock> OperatorType;
>>
>> virtual ~MultiGridTransfer() {}
>>
>> /** \brief Sets up the operator between two P1 spaces
>> *
>> * Does not use any of the Freiburg fem stuff
>> * \param grid The grid
>> * \param cL The coarse grid level
>> * \param fL The fine grid level
>> */
>> template <class GridType>
>> void setup(const GridType& grid, int cL, int fL);
>>
>> /** \brief Restrict a function from the fine onto the coarse grid
>> */
>> virtual void restrict(const DiscFuncType& f, DiscFuncType &t) const;
>>
>> /** \brief Restrict a bitfield from the fine onto the coarse grid
>> */
>> virtual void restrict(const BitField& f, BitField& t) const;
>>
>> /** \brief Prolong a function from the coarse onto the fine grid
>> */
>> virtual void prolong(const DiscFuncType& f, DiscFuncType &t) const;
>>
>> /** \brief Galerkin assemble a coarse stiffness matrix
>> */
>> virtual void galerkinRestrict(const OperatorType& fineMat, OperatorType& coarseMat) const;
>>
>> /** \brief Set Occupation of Galerkin restricted coarse stiffness matrix
>> *
>> * Set occupation of Galerkin restricted coarse matrix. Call this one before
>> * galerkinRestrict to ensure all non-zeroes are present
>> * \param fineMat The fine level matrix
>> * \param coarseMat The coarse level matrix
>> */
>> void galerkinRestrictSetOccupation(const OperatorType& fineMat, OperatorType& coarseMat) const;
>>
>> /** \brief Direct access to the operator matrix, if you absolutely want it! */
>> virtual const OperatorType& getMatrix() const {return matrix_;}
>>
>> protected:
>>
>> OperatorType matrix_;
>>
>> };
>>
>> } // end namespace Dune
>>
>> #include "multigridtransfer.cc"
>>
>> #endif
>
>> #include <dune/common/fvector.hh>
>> #include <dune/istl/matrixindexset.hh>
>> #include <dune/disc/shapefunctions/lagrangeshapefunctions.hh>
>>
>> template<class DiscFuncType>
>> template<class GridType>
>> void Dune::MultiGridTransfer<DiscFuncType>::setup(const GridType& grid, int cL, int fL)
>> {
>> if (fL != cL+1)
>> DUNE_THROW(Exception, "The two function spaces don't belong to consecutive levels!");
>>
>> const int dim = GridType::dimension;
>>
>> const typename GridType::Traits::LevelIndexSet& coarseIndexSet = grid.levelIndexSet(cL);
>> const typename GridType::Traits::LevelIndexSet& fineIndexSet = grid.levelIndexSet(fL);
>>
>> int rows = grid.size(fL, dim);
>> int cols = grid.size(cL, dim);
>>
>> // Make identity matrix
>> MatrixBlock identity(0);
>> for (int i=0; i<blocksize; i++)
>> identity[i][i] = 1;
>>
>> //
>> OperatorType mat(rows, cols, OperatorType::random);
>>
>> mat = 0;
>>
>> typedef typename GridType::template Codim<0>::LevelIterator ElementIterator;
>>
>> ElementIterator cIt = grid.template lbegin<0>(cL);
>> ElementIterator cEndIt = grid.template lend<0>(cL);
>>
>>
>> // ///////////////////////////////////////////
>> // Determine the indices present in the matrix
>> // /////////////////////////////////////////////////
>> MatrixIndexSet indices(rows, cols);
>>
>> for (; cIt != cEndIt; ++cIt) {
>>
>> const LagrangeShapeFunctionSet<double, double, dim>& coarseBaseSet
>> = Dune::LagrangeShapeFunctions<double, double, dim>::general(cIt->geometry().type(), 1);
>> const int numCoarseBaseFct = coarseBaseSet.size();
>>
>> typedef typename GridType::template Codim<0>::Entity EntityType;
>> typedef typename EntityType::HierarchicIterator HierarchicIterator;
>>
>> HierarchicIterator fIt = cIt->hbegin(fL);
>> HierarchicIterator fEndIt = cIt->hend(fL);
>>
>> for (; fIt != fEndIt; ++fIt) {
>>
>> if (fIt->level()==cIt->level())
>> continue;
>>
>> const typename EntityType::Geometry& fGeometryInFather = fIt->geometryInFather();
>>
>> const LagrangeShapeFunctionSet<double, double, dim>& fineBaseSet
>> = Dune::LagrangeShapeFunctions<double, double, dim>::general(fIt->geometry().type(), 1);
>> const int numFineBaseFct = fineBaseSet.size();
>>
>> for (int i=0; i<numCoarseBaseFct; i++) {
>>
>> int globalCoarse = coarseIndexSet.template subIndex<dim>(*cIt,i);
>>
>> for (int j=0; j<numFineBaseFct; j++) {
>>
>> int globalFine = fineIndexSet.template subIndex<dim>(*fIt,j);
>>
>> FieldVector<double, dim> local = fGeometryInFather.global(fineBaseSet[j].position());
>>
>> // Evaluate coarse grid base function
>> double value = coarseBaseSet[i].evaluateFunction(0, local);
>>
>> if (value > 0.001)
>> indices.add(globalFine, globalCoarse);
>>
>> }
>>
>>
>> }
>>
>> }
>>
>> }
>>
>> indices.exportIdx(mat);
>>
>> // /////////////////////////////////////////////
>> // Compute the matrix
>> // /////////////////////////////////////////////
>> cIt = grid.template lbegin<0>(cL);
>> for (; cIt != cEndIt; ++cIt) {
>>
>> const LagrangeShapeFunctionSet<double, double, dim>& coarseBaseSet
>> = Dune::LagrangeShapeFunctions<double, double, dim>::general(cIt->geometry().type(), 1);
>> const int numCoarseBaseFct = coarseBaseSet.size();
>>
>> typedef typename GridType::template Codim<0>::Entity EntityType;
>> typedef typename EntityType::HierarchicIterator HierarchicIterator;
>>
>> HierarchicIterator fIt = cIt->hbegin(fL);
>> HierarchicIterator fEndIt = cIt->hend(fL);
>>
>> for (; fIt != fEndIt; ++fIt) {
>>
>> if (fIt->level()==cIt->level())
>> continue;
>>
>> const typename EntityType::Geometry& fGeometryInFather = fIt->geometryInFather();
>>
>> const LagrangeShapeFunctionSet<double, double, dim>& fineBaseSet
>> = Dune::LagrangeShapeFunctions<double, double, dim>::general(fIt->geometry().type(), 1);
>> const int numFineBaseFct = fineBaseSet.size();
>>
>> for (int i=0; i<numCoarseBaseFct; i++) {
>>
>> int globalCoarse = coarseIndexSet.template subIndex<dim>(*cIt,i);
>>
>> for (int j=0; j<numFineBaseFct; j++) {
>>
>> int globalFine = fineIndexSet.template subIndex<dim>(*fIt,j);
>>
>> // Evaluate coarse grid base function at the location of the fine grid dof
>>
>> // first determine local fine grid dof position
>> FieldVector<double, dim> local = fGeometryInFather.global(fineBaseSet[j].position());
>>
>> // Evaluate coarse grid base function
>> double value = coarseBaseSet[i].evaluateFunction(0, local);
>>
>> // The following conditional is a hack: evaluating the coarse
>> // grid base function will often return 0. However, we don't
>> // want explicit zero entries in our prolongation matrix. Since
>> // the whole code works for P1-elements only anyways, we know
>> // that value can only be 0, 0.5, or 1. Thus testing for nonzero
>> // by testing > 0.001 is safe.
>> if (value > 0.001) {
>> MatrixBlock matValue = identity;
>> matValue *= value;
>> mat[globalFine][globalCoarse] = matValue;
>> }
>>
>> }
>>
>>
>> }
>>
>> }
>>
>> }
>>
>> matrix_ = mat;
>>
>> }
>>
>> // Multiply the vector f from the right to the prolongation matrix
>> template<class DiscFuncType>
>> void Dune::MultiGridTransfer<DiscFuncType>::prolong(const DiscFuncType& f, DiscFuncType& t) const
>> {
>> if (f.size() != matrix_.M())
>> DUNE_THROW(Exception, "Number of entries in the coarse grid vector is not equal "
>> << "to the number of columns of the interpolation matrix!");
>>
>> t.resize(matrix_.N());
>>
>> typedef typename DiscFuncType::Iterator Iterator;
>> typedef typename DiscFuncType::ConstIterator ConstIterator;
>>
>> typedef typename OperatorType::row_type RowType;
>> typedef typename RowType::ConstIterator ColumnIterator;
>>
>> Iterator tIt = t.begin();
>> ConstIterator fIt = f.begin();
>>
>> for(int rowIdx=0; rowIdx<matrix_.N(); rowIdx++) {
>>
>> const RowType& row = matrix_[rowIdx];
>>
>> *tIt = 0.0;
>> ColumnIterator cIt = row.begin();
>> ColumnIterator cEndIt = row.end();
>>
>> for(; cIt!=cEndIt; ++cIt) {
>>
>> cIt->umv(f[cIt.index()], *tIt);
>>
>> }
>>
>> ++tIt;
>> }
>>
>> }
>>
>> // Multiply the vector f from the right to the transpose of the prolongation matrix
>> template<class DiscFuncType>
>> void Dune::MultiGridTransfer<DiscFuncType>::restrict(const DiscFuncType& f, DiscFuncType& t) const
>> {
>> if (f.size() != matrix_.N())
>> DUNE_THROW(Exception, "Fine grid vector has " << f.size() << " entries "
>> << "but the interpolation matrix has " << matrix_.N() << " rows!");
>>
>> t.resize(matrix_.M());
>> t = 0;
>>
>> typedef typename DiscFuncType::Iterator Iterator;
>> typedef typename DiscFuncType::ConstIterator ConstIterator;
>>
>> typedef typename OperatorType::row_type RowType;
>> typedef typename RowType::ConstIterator ColumnIterator;
>>
>> Iterator tIt = t.begin();
>> ConstIterator fIt = f.begin();
>>
>> for (int rowIdx=0; rowIdx<matrix_.N(); rowIdx++) {
>>
>> const RowType& row = matrix_[rowIdx];
>>
>> ColumnIterator cIt = row.begin();
>> ColumnIterator cEndIt = row.end();
>>
>> for(; cIt!=cEndIt; ++cIt) {
>>
>> cIt->umtv(f[rowIdx], t[cIt.index()]);
>>
>> }
>>
>> }
>>
>> }
>>
>>
>> // Multiply the vector f from the right to the transpose of the prolongation matrix
>> template<class DiscFuncType>
>> void Dune::MultiGridTransfer<DiscFuncType>::restrict(const Dune::BitField& f, Dune::BitField& t) const
>> {
>> if (f.size() != (unsigned int)matrix_.N())
>> DUNE_THROW(Exception, "Fine grid bitfield has " << f.size() << " entries "
>> << "but the interpolation matrix has " << matrix_.N() << " rows!");
>>
>> t.resize(matrix_.M());
>> t.unsetAll();
>>
>> typedef typename OperatorType::row_type RowType;
>> typedef typename RowType::ConstIterator ColumnIterator;
>>
>> for (int rowIdx=0; rowIdx<matrix_.N(); rowIdx++) {
>>
>> if (!f[rowIdx])
>> continue;
>>
>> const RowType& row = matrix_[rowIdx];
>>
>> ColumnIterator cIt = row.begin();
>> ColumnIterator cEndIt = row.end();
>>
>> for(; cIt!=cEndIt; ++cIt)
>> t[cIt.index()] = true;
>>
>> }
>>
>> }
>>
>>
>> template<class DiscFuncType>
>> void Dune::MultiGridTransfer<DiscFuncType>::
>> galerkinRestrict(const OperatorType& fineMat, OperatorType& coarseMat) const
>> {
>> // ////////////////////////
>> // Nonsymmetric case
>> // ////////////////////////
>> typedef typename OperatorType::row_type RowType;
>> typedef typename RowType::Iterator ColumnIterator;
>> typedef typename RowType::ConstIterator ConstColumnIterator;
>>
>> // Clear coarse matrix
>> coarseMat = 0;
>>
>> // Loop over all rows of the stiffness matrix
>> for (int v=0; v<fineMat.N(); v++) {
>>
>> const RowType& row = fineMat[v];
>>
>> // Loop over all columns of the stiffness matrix
>> ConstColumnIterator m = row.begin();
>> ConstColumnIterator mEnd = row.end();
>>
>> for (; m!=mEnd; ++m) {
>>
>> int w = m.index();
>>
>> // Loop over all coarse grid vectors iv that have v in their support
>> ConstColumnIterator im = matrix_[v].begin();
>> ConstColumnIterator imEnd = matrix_[v].end();
>> for (; im!=imEnd; ++im) {
>>
>> int iv = im.index();
>>
>> // Loop over all coarse grid vectors jv that have w in their support
>> ConstColumnIterator jm = matrix_[w].begin();
>> ConstColumnIterator jmEnd = matrix_[w].end();
>>
>> for (; jm!=jmEnd; ++jm) {
>>
>> int jv = jm.index();
>>
>> MatrixBlock& cm = coarseMat[iv][jv];
>>
>> // Compute cm = im^T * m * jm
>> for (int i=0; i<blocksize; i++) {
>>
>> for (int j=0; j<blocksize; j++) {
>>
>> for (int k=0; k<blocksize; k++) {
>>
>> for (int l=0; l<blocksize; l++)
>> cm[i][j] += (*im)[k][i] * (*m)[k][l] * (*jm)[l][j];
>>
>> }
>>
>> }
>>
>> }
>>
>> }
>>
>> }
>>
>> }
>>
>> }
>>
>> }
>>
>>
>> // Set Occupation of Galerkin restricted coarse stiffness matrix
>> template<class DiscFuncType>
>> void Dune::MultiGridTransfer<DiscFuncType>::
>> galerkinRestrictSetOccupation(const OperatorType& fineMat, OperatorType& coarseMat) const
>> {
>> // ////////////////////////
>> // Nonsymmetric case
>> // ////////////////////////
>> typedef typename OperatorType::row_type RowType;
>> typedef typename RowType::ConstIterator ConstColumnIterator;
>>
>>
>> // Create index set
>> Dune::MatrixIndexSet indices(matrix_.M(), matrix_.M());
>>
>> // Loop over all rows of the fine matrix
>> for (int v=0; v<fineMat.N(); v++) {
>>
>> const RowType& row = fineMat[v];
>>
>> // Loop over all columns of the fine matrix
>> ConstColumnIterator m = row.begin();
>> ConstColumnIterator mEnd = row.end();
>>
>> for (; m!=mEnd; ++m) {
>>
>> int w = m.index();
>>
>> // Loop over all coarse grid vectors iv that have v in their support
>> ConstColumnIterator im = matrix_[v].begin();
>> ConstColumnIterator imEnd = matrix_[v].end();
>> for (; im!=imEnd; ++im) {
>>
>> int iv = im.index();
>>
>> // Loop over all coarse grid vectors jv that have w in their support
>> ConstColumnIterator jm = matrix_[w].begin();
>> ConstColumnIterator jmEnd = matrix_[w].end();
>>
>> for (; jm!=jmEnd; ++jm)
>> indices.add(iv, jm.index());
>>
>> }
>>
>> }
>>
>> }
>>
>> indices.exportIdx(coarseMat);
>> }
>>
>
>
> --
> Sreejith P. Kuttanikkad
> IWR, University of Heidelberg
> Room:009, Im Neuenheimer Feld-348
> 69120 Heidelberg,Germany.
> Ph :+49-6221-54-5689/4412(Office)
> :+49-(0)17624228904(Mob)
> http://hal.iwr.uni-heidelberg.de/~sreejith/
> -------------------------------------------
>
More information about the Dune
mailing list