From: Benjamin S. Kirk <benkirk@cf...>  20040604 01:46:09

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 ELSheikh 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=(xx_min)/(x_maxx_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 > installationauthoring solution that does it all. Learn more and > evaluate today! http://www.installshield.com/Dev2Dev/0504 > _______________________________________________ > Libmeshdevel mailing list > Libmeshdevel@... > https://lists.sourceforge.net/lists/listinfo/libmeshdevel 
From: Ahmed ELSheikh <elsheiam@mc...>  20040603 15:25:52

Hi, I want to get a patch of all elements and initialize a new submesh with these elements. Say all the neighbors of one element (I didnt use the neighbor iterator because later I will have different types of patchs) So this is what I did, //// const_active_local_elem_iterator elem_it (mesh.elements_begin()); const const_active_local_elem_iterator elem_end(mesh.elements_end()); for (; elem_it != elem_end; ++elem_it) { const Elem* input_el= *elem_it; Mesh patch_mesh (2); // a vector to store the patch elements std::vector<const Elem *> elements_patch; elements_patch.reserve(input_el>n_neighbors()+1); // insert the element itself elements_patch.push_back(input_el); for (unsigned int n_e=0; n_e<input_el>n_neighbors(); n_e++) if (input_el>neighbor(n_e) != NULL) { Elem* ef = input_el>neighbor(n_e); elements_patch.push_back(ef); } // end neighbors // create iterators to pass it to create_submesh const_elem_iterator it(std::make_pair( elements_patch.begin(), elements_patch.end())); const const_elem_iterator it_end(std::make_pair( ((std::vector<const Elem*>*) &elements_patch)>end(), ((std::vector<const Elem*>*) &elements_patch)>end())); mesh.create_submesh (patch_mesh, it, it_end); patch_mesh.print_info(); } // end of elem_it // 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*, std::allocator<const Elem*> > >&) /usr/local/include/c++/3.3.3/bits/stl_iterator.h:593: error: __gnu_cxx::__normal_iterator<_Iterator, _Container>::__normal_iterator(const _Iterator&) [with _Iterator = const Elem* const*, _Container = std::vector<const Elem*, std::allocator<const Elem*> >] /usr/local/include/c++/3.3.3/bits/stl_iterator.h:590: error: __gnu_cxx::__normal_iterator<_Iterator, _Container>::__normal_iterator() [with _Iterator = const Elem* const*, _Container = std::vector<const Elem*, std::allocator<const Elem*> >] /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*, std::allocator<const Elem*> > >&) /usr/local/include/c++/3.3.3/bits/stl_iterator.h:593: error: __gnu_cxx::__normal_iterator<_Iterator, _Container>::__normal_iterator(const _Iterator&) [with _Iterator = const Elem* const*, _Container = std::vector<const Elem*, std::allocator<const Elem*> >] /usr/local/include/c++/3.3.3/bits/stl_iterator.h:590: error: __gnu_cxx::__normal_iterator<_Iterator, _Container>::__normal_iterator() [with _Iterator = const Elem* const*, _Container = std::vector<const Elem*, std::allocator<const Elem*> >] 
From: John Peterson <peterson@cf...>  20040603 17:17:51

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 differentlynamed function. I hope that helps you. John Ahmed ELSheikh 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*, 
From: Ahmed ELSheikh <elsheiam@mc...>  20040603 17:57:09

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. ************************************************** 
From: John Peterson <peterson@cf...>  20040603 18:08:01

Ahmed ELSheikh 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 ZienkiewiczZhu superconvergent patch recovery error indicator which computes an improved approximation to the flux grad(u) on patches of userdefinable 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 
From: Ahmed ELSheikh <elsheiam@mc...>  20040603 18:35:53

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=(xx_min)/(x_maxx_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 
From: John Peterson <peterson@cf...>  20040603 19:52:21

Ahmed ELSheikh 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=(xx_min)/(x_maxx_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 residualbased error indicators, but by definition they are problemspecific and not generally appropriate for inclusion in the library. The best approach (as you mention below) is to reuse 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 interelement 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 
From: Benjamin S. Kirk <benkirk@cf...>  20040604 01:46:09

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 ELSheikh 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=(xx_min)/(x_maxx_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 > installationauthoring solution that does it all. Learn more and > evaluate today! http://www.installshield.com/Dev2Dev/0504 > _______________________________________________ > Libmeshdevel mailing list > Libmeshdevel@... > https://lists.sourceforge.net/lists/listinfo/libmeshdevel 
From: Ahmed ELSheikh <elsheiam@mc...>  20040604 04:08:12
Attachments:
svd_sent_to_ben.h

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 ELSheikh 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=(xx_min)/(x_maxx_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 > > installationauthoring solution that does it all. Learn more and > > evaluate today! http://www.installshield.com/Dev2Dev/0504 > > _______________________________________________ > > Libmeshdevel mailing list > > Libmeshdevel@... > > https://lists.sourceforge.net/lists/listinfo/libmeshdevel 
From: John Peterson <peterson@cf...>  20040604 14:44:17

Ahmed ELSheikh 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 