libMesh: A C++ Finite Element Library

Hi,

Neat idea.  Does your submesh have a hole in it, or do you
insert the current element as well?

Try the following line:

const_active_local_elem_iterator       pelem_it
(patch_mesh.const_elements_begin());
            *****

const const_active_local_elem_iterator
pelem_end (patch_mesh.const_elements_end());
                      *****


Since patch_mesh is not a const object, but you want const_iterators
out of it, you have to use the differently-named function.
I hope that helps you.

-John






Ahmed EL-Sheikh writes:
 > 
 > So far everything works fine and the patch_mesh.print_info(); shows a
 > correct size of the patch.
 > 
 > When I try to extract infromation from the patch_mesh I used
 > 
 > const_active_local_elem_iterator       pelem_it
 > (patch_mesh.elements_begin());
 > const const_active_local_elem_iterator
 > pelem_end(patch_mesh.elements_end());
 > 
 > for (; pelem_it != pelem_end; ++pelem_it)
 > {
 > const Elem* e = *pelem_it;
 > e->write_tecplot_connectivity(std::cout);
 > for (unsigned int n=0; n<e->n_nodes(); n++)   e->point(n).print();
 > 
 > }
 > 
 > But I got an error related to the iterators. I'll attach it at the end
 > of
 > this email.
 > 
 > Any advice of what I am missing.
 > 
 > Ahmed
 > 
 > 
 > ///// ------------------- compilation errors
 > /usr/local/include/c++/3.3.3/bits/stl_pair.h: In constructor
 > `std::pair<_T1,
 >    _T2>::pair(const std::pair<_U1, _U2>&) [with _U1 =
 >    __gnu_cxx::__normal_iterator<Elem**, std::vector<Elem*,
 >    std::allocator<Elem*> > >, _U2 = __gnu_cxx::__normal_iterator<Elem**,
 > 
 >    std::vector<Elem*, std::allocator<Elem*> > >, _T1 =
 >    __gnu_cxx::__normal_iterator<const Elem* const*, std::vector<const
 > Elem*,
 >    std::allocator<const Elem*> > >, _T2 =
 > __gnu_cxx::__normal_iterator<const
 >    Elem* const*, std::vector<const Elem*, std::allocator<const Elem*> >
 > >]':
 > /work/elsheiam/libmesh/include/mesh/elem_iterators.h:851:   instantiated
 > 
 > from here
 > /usr/local/include/c++/3.3.3/bits/stl_pair.h:88: error: no matching
 > function
 >    for call to `__gnu_cxx::__normal_iterator<const Elem* const*,
 >    std::vector<const Elem*, std::allocator<const Elem*> >
 > >::__normal_iterator(
 >    const __gnu_cxx::__normal_iterator<Elem**, std::vector<Elem*,
 >    std::allocator<Elem*> > >&)'
 > /usr/local/include/c++/3.3.3/bits/stl_iterator.h:580: error: candidates
 > are:
 >    __gnu_cxx::__normal_iterator<const Elem* const*, std::vector<const
 > Elem*,
 >    std::allocator<const Elem*> > >::__normal_iterator(const
 >    __gnu_cxx::__normal_iterator<const Elem* const*, std::vector<const
 > Elem*,

Thanks John so much for this fast reply. It works :)

Actually, I am investigating another error estimation
technique based on solving local patches.

Is anybody interested or did similar work before?

> Neat idea.  Does your submesh have a hole in it, or do you
> insert the current element as well?

I have already inserted the current element before iterating on the neighbors.
Anyway, I will check that again.

Regards,
Ahmed

--
"No great thing is created suddenly." Epictetus
**************************************************
  Ahmed ElSheikh <elsheiam_AT_mcmaster.ca>
  PhD Candidate, Department of Civil Engineering
  McMaster University
  1280 Main Street West
  Hamilton, Ontario,  L8S 4L7, Canada.
**************************************************

Ahmed EL-Sheikh writes:
 > Thanks John so much for this fast reply. It works :)
 > 
 > Actually, I am investigating another error estimation
 > technique based on solving local patches.
 > 
 > Is anybody interested or did similar work before?

Actually, yes, we have just been discussing a generalization of
the Zienkiewicz-Zhu superconvergent patch recovery error indicator
which computes an improved approximation to the flux grad(u)
on patches of user-definable size.  

Note that your technique will only add *face* neighbors of the
current element to the patch, i.e. you will not in general include
the full support of each of the nodal basis functions in your
patch.  The latest CVS should have the start of the patch recovery error 
code in it, if you want to take a look.

 > > Neat idea.  Does your submesh have a hole in it, or do you
 > > insert the current element as well?
 > 
 > I have already inserted the current element before iterating on the neighbors.
 > Anyway, I will check that again.

-John

That is great.. I have already developed a similar code :-(
I started with the same code but I added some modifications

The first is the need to use a local coordinate system for the
patch to avoid the different weighting of   x and y.
So I did something like this:

// determine the x_min, x_max, y_min, y_max , z_min, z_max for all the points on the
patch
// Scaled x and y by z by
// x=(x-x_min)/(x_max-x_min);


The other thing is using Singular value decomposition to solve.
I adopted a code from TNT for that.

Do you have any plans for different error estimators ....


John Peterson wrote:

> Note that your technique will only add *face* neighbors of the
> current element to the patch, i.e. you will not in general include
> the full support of each of the nodal basis functions in your
> patch.

I new that .. It is a bit problem here,
What I want is a patch ... with the same system of equation
as well as the current solution. This patch may use the same assemble function
but on the residuals or attach its own assembling function.


One choice is use these faces I build a map from the patch mesh to
the original mesh ... then initialize a different equation system for the
patch .. obtain what is needed from the old mesh when needed.
Attach another assembling function for the error estimation with
different boundary conditions. Then instead of solving it using the method
solve use SVD.

Is that make sense !!

Ahmed

Ahmed EL-Sheikh writes:
 > That is great.. I have already developed a similar code :-(
 > I started with the same code but I added some modifications
 > 
 > The first is the need to use a local coordinate system for the
 > patch to avoid the different weighting of   x and y.
 > So I did something like this:
 > 
 > // determine the x_min, x_max, y_min, y_max , z_min, z_max for all the points on the
 > patch
 > // Scaled x and y by z by
 > // x=(x-x_min)/(x_max-x_min);
 > 
 > 
 > The other thing is using Singular value decomposition to solve.
 > I adopted a code from TNT for that.
 > 
 > Do you have any plans for different error estimators ....

I have been thinking briefly about residual-based error indicators, but
by definition they are problem-specific and not generally appropriate
for inclusion in the library.  The best approach (as you mention below)
is to re-use the same assemble function already provided by the user
to the EquationSystems, but do not assemble any global matrix, just
compute the residual on each element.  That might take some modification
to the way we currently code the assembly routines, but it could be done.

A problem with the above approach is that the edge residuals,
i.e. the jump in the normal derivatives across inter-element boundaries,
are sometimes the largest contribution to an element residual, and they
are not computed during most matrix assembly routines (for continuous
Galerkin).  So providing a generic

class ResidualBasedErrorEstimator {};

with the method

attach_residual_function(...);

would probably be what would eventually make it into the library.  

 > One choice is use these faces I build a map from the patch mesh to
 > the original mesh ... then initialize a different equation system for the
 > patch .. obtain what is needed from the old mesh when needed.
 > Attach another assembling function for the error estimation with
 > different boundary conditions. Then instead of solving it using the method
 > solve use SVD.

Hmm...are the systems of equations that you generate singular or nearly singular?
We assumed that it would be possible to compute an LU decomposition to solve
for the coefficients, but maybe that won't work....?

-John

What a relevant thread!  Funny, as John mentioned, we started 
implementing this yesterday.

What approximation order are you using for the reconstructed gradient? 
Certainly, when the patch size is close to the number of unknowns in the 
reconstruction (or less!) you can get a horribly conditioned (or 
singular) system that needs SVD.  We have some numerical experiments 
that suggest LU is OK for linear and quadratic reconstructions when the 
patches are of "reasonable" size...

Eventually I'd like the library to look for LAPACK in ./configure and 
then use it internal to the DenseMatrix<> and DenseVector<> class.  In 
the meantime, if you'd like to "share" an SVD implementation we'd be 
happy to include it!

-Ben

Ahmed EL-Sheikh wrote:
> That is great.. I have already developed a similar code :-(
> I started with the same code but I added some modifications
> 
> The first is the need to use a local coordinate system for the
> patch to avoid the different weighting of   x and y.
> So I did something like this:
> 
> // determine the x_min, x_max, y_min, y_max , z_min, z_max for all the points on the
> patch
> // Scaled x and y by z by
> // x=(x-x_min)/(x_max-x_min);
> 
> 
> The other thing is using Singular value decomposition to solve.
> I adopted a code from TNT for that.
> 
> Do you have any plans for different error estimators ....
> 
> 
> John Peterson wrote:
> 
> 
>>Note that your technique will only add *face* neighbors of the
>>current element to the patch, i.e. you will not in general include
>>the full support of each of the nodal basis functions in your
>>patch.
> 
> 
> I new that .. It is a bit problem here,
> What I want is a patch ... with the same system of equation
> as well as the current solution. This patch may use the same assemble function
> but on the residuals or attach its own assembling function.
> 
> 
> One choice is use these faces I build a map from the patch mesh to
> the original mesh ... then initialize a different equation system for the
> patch .. obtain what is needed from the old mesh when needed.
> Attach another assembling function for the error estimation with
> different boundary conditions. Then instead of solving it using the method
> solve use SVD.
> 
> Is that make sense !!
> 
> Ahmed
> 
> 
> 
> -------------------------------------------------------
> This SF.Net email is sponsored by the new InstallShield X.
>>From Windows to Linux, servers to mobile, InstallShield X is the one
> installation-authoring solution that does it all. Learn more and
> evaluate today! http://www.installshield.com/Dev2Dev/0504
> _______________________________________________
> Libmesh-devel mailing list
> Libmesh-devel@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-devel

The problem with the number of unknowns appears for the elements on the boundary.. linear
chain of quad elements is simple example. As you said for linear and quadratic
reconstructions it is some how ok but if the initial mesh is bad there might be a problem.

> Eventually I'd like the library to look for LAPACK in ./configure and
> then use it internal to the DenseMatrix<> and DenseVector<> class.  In
> the meantime, if you'd like to "share" an SVD implementation we'd be
> happy to include it!

I would recommend to adopt an SVD implementation like the LU one so if you run it in
parallel it works too. I will attach the code which I have it is the TNT implementation of
the SVD. I modified it to work with the DenseMatrix<> and DenseVector<>.
I call it from the error estimator like

//----------------------------
SVD<double> M_svd(M);
M_svd.solve(bx, solx);
M_svd.solve(by, soly);
//----------------------------

BUT this code may loop for ever. I am checking how to put a stopping criteria if there is
a convergence problem.

Ahmed


>
> -Ben
>
> Ahmed EL-Sheikh wrote:
> > That is great.. I have already developed a similar code :-(
> > I started with the same code but I added some modifications
> >
> > The first is the need to use a local coordinate system for the
> > patch to avoid the different weighting of   x and y.
> > So I did something like this:
> >
> > // determine the x_min, x_max, y_min, y_max , z_min, z_max for all the points on the
> > patch
> > // Scaled x and y by z by
> > // x=(x-x_min)/(x_max-x_min);
> >
> >
> > The other thing is using Singular value decomposition to solve.
> > I adopted a code from TNT for that.
> >
> > Do you have any plans for different error estimators ....
> >
> >
> > John Peterson wrote:
> >
> >
> >>Note that your technique will only add *face* neighbors of the
> >>current element to the patch, i.e. you will not in general include
> >>the full support of each of the nodal basis functions in your
> >>patch.
> >
> >
> > I new that .. It is a bit problem here,
> > What I want is a patch ... with the same system of equation
> > as well as the current solution. This patch may use the same assemble function
> > but on the residuals or attach its own assembling function.
> >
> >
> > One choice is use these faces I build a map from the patch mesh to
> > the original mesh ... then initialize a different equation system for the
> > patch .. obtain what is needed from the old mesh when needed.
> > Attach another assembling function for the error estimation with
> > different boundary conditions. Then instead of solving it using the method
> > solve use SVD.
> >
> > Is that make sense !!
> >
> > Ahmed
> >
> >
> >
> > -------------------------------------------------------
> > This SF.Net email is sponsored by the new InstallShield X.
> >>From Windows to Linux, servers to mobile, InstallShield X is the one
> > installation-authoring solution that does it all. Learn more and
> > evaluate today! http://www.installshield.com/Dev2Dev/0504
> > _______________________________________________
> > Libmesh-devel mailing list
> > Libmesh-devel@...
> > https://lists.sourceforge.net/lists/listinfo/libmesh-devel

Ahmed EL-Sheikh writes:
 > The problem with the number of unknowns appears for the elements on the boundary.. linear
 > chain of quad elements is simple example. As you said for linear and quadratic
 > reconstructions it is some how ok but if the initial mesh is bad there might be a problem.
 > 
 > > Eventually I'd like the library to look for LAPACK in ./configure and
 > > then use it internal to the DenseMatrix<> and DenseVector<> class.  In
 > > the meantime, if you'd like to "share" an SVD implementation we'd be
 > > happy to include it!
 > 
 > I would recommend to adopt an SVD implementation like the LU one so if you run it in
 > parallel it works too. I will attach the code which I have it is the TNT implementation of
 > the SVD. I modified it to work with the DenseMatrix<> and DenseVector<>.
 > I call it from the error estimator like

FYI: Here is the Copyright notice from the page if anyone is interested.
     It's public domain software, so I guess there should be no problem copying/distributing
     it under LGPL.




Copyright Notice This software is a cooperative product of The
MathWorks and the National Institute of Standards and Technology
(NIST) which has been released to the public domain. Neither The
MathWorks nor NIST assumes any responsibility whatsoever for its use
by other parties, and makes no guarantees, expressed or implied, about
its quality, reliability, or any other characteristic.


-John