From: Min Xu <minxu@sc...>  20050912 21:41:56

Hi, I am a new libmesh user. I try to solve an optimization problem for the absorption alpha, governed by the diffusion equation: dpsi/dt = laplace(psi)  alpha*psi, such that the psi on the boundary matches the measurement. The absorption alpha is space dependent and need to be determined through optimization. The initial guess for alpha is uniform. I would like to use a coarse mesh for this parameter alpha initially and to refine this mesh gradually in the optimization process at locations where alpha changes abruptly. The solution psi will be solved on a different mesh which may start from the coefficient mesh and evolve independently. I've read the available documents and am convinced this idea is possible to implement by libmesh. I am unclear though what is the way to implement it. Can I just define a dummy EquationSystems to hold this parameter alpha and treat it as "solution", so all the operations on solution (mesh refinement and projection) can be used for alpha and then mapped to the diffusion equationSystems? Is there better ways to achieve it? Thanks and best regards, Min Xu 
From: John Peterson <peterson@cf...>  20050913 01:22:32

Min Xu writes: > Hi, > > I am a new libmesh user. I try to solve an optimization problem for the > absorption alpha, governed by the diffusion equation: > dpsi/dt = laplace(psi)  alpha*psi, > such that the psi on the boundary matches the measurement. The absorption > alpha is space dependent and need to be determined through optimization. > > The initial guess for alpha is uniform. I would like to use a coarse mesh for > this parameter alpha initially and to refine this mesh gradually in the > optimization process at locations where alpha changes abruptly. > > The solution psi will be solved on a different mesh which may start from the > coefficient mesh and evolve independently. I wonder will you really need two separate meshes? It seems that psi will vary a lot wherever alpha varies a lot  they are coupled. The one mesh case I think is simpler than having multiple meshes. Then alpha is similar to a "lagged" coefficient which gives you a new value of psi that in turn feeds back in to determine the next alpha. John 
From: Min Xu <minxu@sc...>  20050913 02:41:10

The solution psi will also vary significantly nearby the source or detector whereas the absorption alpha may just be uniform on these places. It is hence beneficiary to use a more coarse mesh for alpha than that for psi. I am now leaning toward MeshData to store alpha. I read the mailing list archive that MeshData did not support MeshRefinement yet. Have this situation changed? Is it possible to use MeshRefinement to refine the MeshData in the same way as the solution of an equation system is refined? How? Thanks! Min On Monday 12 September 2005 21:22, John Peterson wrote: > Min Xu writes: > > Hi, > > > > I am a new libmesh user. I try to solve an optimization problem for the > > absorption alpha, governed by the diffusion equation: > > dpsi/dt = laplace(psi)  alpha*psi, > > such that the psi on the boundary matches the measurement. The > > absorption alpha is space dependent and need to be determined through > > optimization. > > > > The initial guess for alpha is uniform. I would like to use a coarse > > mesh for this parameter alpha initially and to refine this mesh > > gradually in the optimization process at locations where alpha changes > > abruptly. > > > > The solution psi will be solved on a different mesh which may start from > > the coefficient mesh and evolve independently. > > I wonder will you really need two separate meshes? It seems that psi will > vary a lot wherever alpha varies a lot  they are coupled. The one mesh > case I think is simpler than having multiple meshes. Then alpha is similar > to a "lagged" coefficient which gives you a new value of psi that in turn > feeds back in to determine the next alpha. > > John 
From: John Peterson <peterson@cf...>  20050913 03:36:53

Min Xu writes: > The solution psi will also vary significantly nearby the source or detector > whereas the absorption alpha may just be uniform on these places. It is hence > beneficiary to use a more coarse mesh for alpha than that for psi. That may be true, but you will save yourself a lot of trouble and get started faster using a single mesh. > I am now leaning toward MeshData to store alpha. I read the mailing list > archive that MeshData did not support MeshRefinement yet. Have this situation > changed? Is it possible to use MeshRefinement to refine the MeshData in the > same way as the solution of an equation system is refined? How? I wish I knew more about MeshData, unfortunately I've never used it. John 
From: Roy Stogner <roystgnr@ic...>  20050913 04:36:46

On Mon, 12 Sep 2005, Min Xu wrote: > Can I just define a dummy EquationSystems to hold this parameter alpha and > treat it as "solution", so all the operations on solution (mesh refinement > and projection) can be used for alpha and then mapped to the diffusion > equationSystems? Is there better ways to achieve it? This will work, but you'll have to have a separate mesh for each of the two EquationSystems objects and you'll have to have some efficient way of mapping between cells of one mesh and cells of the other. This map could be two vectors, alpha_to_psi and psi_to_alpha, which each are indexed by element ids from one mesh and return as data element ids from the other. Unfortunately these ids are renumbered every time the mesh is changed. So, after every EquationSystems::reinit() you'd need to regenerate your map vectors from scratch, which could be an expensive process unless you have a very coarse initial grid. Ideally the MeshData class could be used to preserve enough data through refinement to make this regeneration O(N), but like John I don't use MeshData myself.  Roy Stogner 
From: Min Xu <minxu@sc...>  20050913 17:17:04

I am thinking to attach two system to one equationsystem like: EquationSystems equation_systems (mesh); // The real system for diffusion equation TransientLinearImplicitSystem & system = equation_systems.add_system<TransientLinearImplicitSystem>("Diffusion"); system.add_variable ("u", FIRST); // The dummy system for the coefficient alpha LinearImplicitSystem & dummy = equation_systems.add_system<LinearImplicitSystem>("Dummy"); system.add_variable ("alpha", FIRST); The procedure follows: (1) assign a coarse mesh and initialize the systems, mesh_psi = mesh_alpha = mesh0. keep a copy of the mesh as mesh_coarse. (2) solve the diffusion system and allow mesh_psi to be refined in its solution. obtain the update dalpha defined on mesh_psi_refined. (3) project dalpha to mesh_coarse which is the parent mesh of mesh_psi_refined. (4) Update alpha to alpha+dalpha, and refine the mesh for the dummy system to mesh_new_coarse. (5) Set mesh_psi=mesh_alpha=mesh_coarse=mesh_new_coarse and return to (2). My questions are: On Tuesday 13 September 2005 00:36, Roy Stogner wrote: > On Mon, 12 Sep 2005, Min Xu wrote: > > Can I just define a dummy EquationSystems to hold this parameter alpha > > and treat it as "solution", so all the operations on solution (mesh > > refinement and projection) can be used for alpha and then mapped to the > > diffusion equationSystems? Is there better ways to achieve it? > > This will work, but you'll have to have a separate mesh for each of > the two EquationSystems objects and you'll have to have some efficient > way of mapping between cells of one mesh and cells of the other. This > map could be two vectors, alpha_to_psi and psi_to_alpha, which each > are indexed by element ids from one mesh and return as data element > ids from the other. > > Unfortunately these ids are renumbered every time the mesh is changed. > So, after every EquationSystems::reinit() you'd need to regenerate > your map vectors from scratch, which could be an expensive process > unless you have a very coarse initial grid. Ideally the MeshData > class could be used to preserve enough data through refinement to make > this regeneration O(N), but like John I don't use MeshData myself. >  (1) If the two meshes are strictly that one is a refined version of the another (maybe refined multiple times), is there a relative cheap way to use alpha defined on the coarser mesh in the assemble routine for the equation system? Any example for coefficients defined on a mesh? (2) Can project_vector work in a reverse direction, i.e., map a vector defined on a finer mesh to a coarser mesh? If not, any suggestion to implement such a projection? Thanks! Min 
From: John Peterson <peterson@cf...>  20050913 17:43:24

Min Xu writes: > > (1) If the two meshes are strictly that one is a refined version of the > another (maybe refined multiple times), is there a relative cheap way to use > alpha defined on the coarser mesh in the assemble routine for the equation > system? Any example for coefficients defined on a mesh? Not that I know of. In general two separate mesh objects MeshA and MeshB are entirely unrelated, and there's no way to tell them that they are related. If MeshA is a coarse submesh of MeshB, then theoretically it would be possible to coarsen the elements of MeshB in just the right way to obtain a Mesh identical to MeshA. (You would have to reproject the solution at each coarsening step, this wouldn't be "relatively cheap".) Seriously, just use a single mesh for all the variables :) > (2) Can project_vector work in a reverse direction, i.e., map a vector defined > on a finer mesh to a coarser mesh? If not, any suggestion to implement such a > projection? I don't think it can...as far as I know it only projects old vectors onto the current mesh. John 
From: Roy Stogner <roystgnr@ic...>  20050913 17:45:51

On Tue, 13 Sep 2005, Min Xu wrote: > I am thinking to attach two system to one equationsystem like: > > EquationSystems equation_systems (mesh); > > // The real system for diffusion equation > TransientLinearImplicitSystem & system = > equation_systems.add_system<TransientLinearImplicitSystem>("Diffusion"); > system.add_variable ("u", FIRST); > > // The dummy system for the coefficient alpha > LinearImplicitSystem & dummy = > equation_systems.add_system<LinearImplicitSystem>("Dummy"); > system.add_variable ("alpha", FIRST); This is one way to do things if you want both alpha and u to be defined on the same mesh, however unless you're going to be doing implicit solves for alpha you might want to make dummy an explicit system and avoid allocating an extra matrix. > The procedure follows: > (1) assign a coarse mesh and initialize the systems, mesh_psi = > mesh_alpha = mesh0. keep a copy of the mesh as mesh_coarse. You can keep a copy of the mesh, but alpha will always be defined on the mesh attached to equation_systems. If you want two different meshes for alpha and psi, there's currently no way to do so without two separate EquationSystems objects. > (1) If the two meshes are strictly that one is a refined version of the > another (maybe refined multiple times), is there a relative cheap way to use > alpha defined on the coarser mesh in the assemble routine for the equation > system? Any example for coefficients defined on a mesh? You're right that the problem is easier if one mesh is a superset of the other: if your alpha mesh is a coarser version of your psi mesh then a normal quadrature rule on psi will still be accurate when integrating terms containing alpha. However, the trick from a coding standpoint remains the same: in your assemble routine, you'll iterate through elements on the psi mesh. You'll need to integrate against data from the alpha mesh. This means you'll need to be able first to take a mesh_psi element and find the corresponding mesh_alpha element (via some map like I mentioned earlier) and second to take a mesh_psi master element point and find the corresponding mesh_alpha master element point (via inverse_map(), if you want to start with inefficient but correct and quicktoprogram code). > (2) Can project_vector work in a reverse direction, i.e., map a > vector defined on a finer mesh to a coarser mesh? If not, any > suggestion to implement such a projection? Currently the meshtomesh project_vector only works in the context of a single mesh whose elements are being refined and/or coarsened one level at a time, and so it's only meant to be called internally by the EquationSystems object. You could project all your EquationSystems vectors from a finer mesh to a much coarser mesh by repeatedly coarsening and calling EquationSystems::reinit(), but this is a one way process  you wouldn't be able to save fine mesh data. My suggestion would be to do such a projection with the same two EquationSystems objects you'd want to use for the rest of your problem. Either make alpha an implicit system and do an L2 projection, or simply do interpolation. Doing a projection would be the harder to code of the two choices, but would give you code portable to nonLagrange elements. Doing interpolation with Lagrange elements should be easy: iterate through all your mesh_alpha elements, find the corresponding mesh_phi elements, and copy over values. The mesh_phi elements will be inactive "ancestors", but they'll have the nodal data you want.  Roy Stogner 
From: Min Xu <minxu@sc...>  20050914 15:26:10

On Tuesday 13 September 2005 13:45, Roy Stogner wrote: > On Tue, 13 Sep 2005, Min Xu wrote: > > I am thinking to attach two system to one equationsystem like: > > > > EquationSystems equation_systems (mesh); > > > > // The real system for diffusion equation > > TransientLinearImplicitSystem & system = > > > > equation_systems.add_system<TransientLinearImplicitSystem>("Diffusion"); > > system.add_variable ("u", FIRST); > > > > // The dummy system for the coefficient alpha > > LinearImplicitSystem & dummy = > > equation_systems.add_system<LinearImplicitSystem>("Dummy"); > > system.add_variable ("alpha", FIRST); > > This is one way to do things if you want both alpha and u to be > defined on the same mesh, however unless you're going to be doing > implicit solves for alpha you might want to make dummy an explicit > system and avoid allocating an extra matrix. > > > The procedure follows: > > (1) assign a coarse mesh and initialize the systems, mesh_psi = > > mesh_alpha = mesh0. keep a copy of the mesh as mesh_coarse. > > You can keep a copy of the mesh, but alpha will always be defined on > the mesh attached to equation_systems. If you want two different > meshes for alpha and psi, there's currently no way to do so without > two separate EquationSystems objects. > You are right. I could imagine though it may be desirable to allow the multiple systems defined within one equationSystem to use their individual meshes. For example, L1 psi1 = f (psi2) L2 psi2 = g(psi1) one can use different meshes for psi1 and psi2 as seen most effective. > > (1) If the two meshes are strictly that one is a refined version of the > > another (maybe refined multiple times), is there a relative cheap way to > > use alpha defined on the coarser mesh in the assemble routine for the > > equation system? Any example for coefficients defined on a mesh? > > You're right that the problem is easier if one mesh is a superset of > the other: if your alpha mesh is a coarser version of your psi mesh > then a normal quadrature rule on psi will still be accurate when > integrating terms containing alpha. > > However, the trick from a coding standpoint remains the same: in your > assemble routine, you'll iterate through elements on the psi mesh. > You'll need to integrate against data from the alpha mesh. This means > you'll need to be able first to take a mesh_psi element and find the > corresponding mesh_alpha element (via some map like I mentioned > earlier) and second to take a mesh_psi master element point and find > the corresponding mesh_alpha master element point (via inverse_map(), > if you want to start with inefficient but correct and quicktoprogram > code). The corresponding element number for the same spatial location may be different between the two meshes due to renumbering as you have pointed out in the earlier email. To implement such a mapping, how should I proceed? What is the index used for the solution vector? > > > (2) Can project_vector work in a reverse direction, i.e., map a > > vector defined on a finer mesh to a coarser mesh? If not, any > > suggestion to implement such a projection? > > Currently the meshtomesh project_vector only works in the context of > a single mesh whose elements are being refined and/or coarsened one > level at a time, and so it's only meant to be called internally by the > EquationSystems object. You could project all your EquationSystems > vectors from a finer mesh to a much coarser mesh by repeatedly > coarsening and calling EquationSystems::reinit(), but this is a one > way process  you wouldn't be able to save fine mesh data. > > My suggestion would be to do such a projection with the same two > EquationSystems objects you'd want to use for the rest of your > problem. Either make alpha an implicit system and do an L2 > projection, or simply do interpolation. Doing a projection would be > the harder to code of the two choices, but would give you code > portable to nonLagrange elements. Doing interpolation with Lagrange > elements should be easy: iterate through all your mesh_alpha elements, > find the corresponding mesh_phi elements, and copy over values. The > mesh_phi elements will be inactive "ancestors", but they'll have the > nodal data you want. The required functionality conceptually belongs to MeshData. However, I do not find any way to use MeshData for this purpose such that data defined on a mesh can be casted into a new MeshData on its child or parent mesh by projection or restriction. Anyone knows a way? Thanks! Min 
From: Roy Stogner <roystgnr@ic...>  20050914 16:07:28

On Wed, 14 Sep 2005, Min Xu wrote: > You are right. I could imagine though it may be desirable to allow the > multiple systems defined within one equationSystem to use their individual > meshes. For example, > L1 psi1 = f (psi2) > L2 psi2 = g(psi1) > one can use different meshes for psi1 and psi2 as seen most effective. I wrote code (using deal.II) to try and do this once, and I agree that it's desirable in many cases, but it was a failure for me. You have to be careful doing this with mixed methods; everyone knows about the element stability restrictions on velocity/pressure elements, but it's not as obvious that even systems like streamfunction/vorticity can become over or underdetermined if you try to refine the solution components independently. > The corresponding element number for the same spatial location may be > different between the two meshes due to renumbering as you have pointed out > in the earlier email. To implement such a mapping, how should I proceed? Start by setting up the mapping on the coarsest elements, for which there should be a onetoone correspondence. I don't know of any way to do this in general that isn't O(n^2), but if you're willing to risk your code being broken by future libmesh revisions, you can make the assumption that coarse element ids don't change, and so "map[i] == i" for 0 < i < num_coarse_elements will remain true no matter how you refine the two meshes. Then, for each coarse element, do a depthfirst search down both element's refinement trees. If elemA and elemB map to each other, then so will elemA>child(i) and elemB>child(i). If elemB has children but elemA does not, then all of elemB's children (and grandchildren, etc) will map to elem(i). > What is the index used for the solution vector? The indexing will depend on the system. Each system will have its own separate dof_map which you'll need to convert the local indices into global indices.  Roy Stogner 
From: Min Xu <minxu@sc...>  20050915 11:51:55

It turns out the functionality requested _has_ been implemented in MeshFunction. A while ago, Michael Schindler sent in a patch for MeshFunction and other minor code corrections to this list. Has the patch been merged into the official libmesh? Min On Wednesday 14 September 2005 12:07, Roy Stogner wrote: > On Wed, 14 Sep 2005, Min Xu wrote: > > You are right. I could imagine though it may be desirable to allow the > > multiple systems defined within one equationSystem to use their > > individual meshes. For example, > > L1 psi1 = f (psi2) > > L2 psi2 = g(psi1) > > one can use different meshes for psi1 and psi2 as seen most effective. > > I wrote code (using deal.II) to try and do this once, and I agree that > it's desirable in many cases, but it was a failure for me. You have > to be careful doing this with mixed methods; everyone knows about the > element stability restrictions on velocity/pressure elements, but it's > not as obvious that even systems like streamfunction/vorticity can > become over or underdetermined if you try to refine the solution > components independently. > > > The corresponding element number for the same spatial location may be > > different between the two meshes due to renumbering as you have pointed > > out in the earlier email. To implement such a mapping, how should I > > proceed? > > Start by setting up the mapping on the coarsest elements, for which > there should be a onetoone correspondence. I don't know of any way > to do this in general that isn't O(n^2), but if you're willing to risk > your code being broken by future libmesh revisions, you can make the > assumption that coarse element ids don't change, and so "map[i] == i" > for 0 < i < num_coarse_elements will remain true no matter how you > refine the two meshes. > > Then, for each coarse element, do a depthfirst search down both > element's refinement trees. If elemA and elemB map to each other, > then so will elemA>child(i) and elemB>child(i). If elemB has > children but elemA does not, then all of elemB's children (and > grandchildren, etc) will map to elem(i). > > > What is the index used for the solution vector? > > The indexing will depend on the system. Each system will have its > own separate dof_map which you'll need to convert the local indices > into global indices. >  > Roy Stogner > > >  > SF.Net email is sponsored by: > Tame your development challenges with Apache's Geronimo App Server. > Download it for free  and be entered to win a 42" plasma tv or your very > own Sony(tm)PSP. Click here to play: http://sourceforge.net/geronimo.php > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers 
From: Roy Stogner <roystgnr@ic...>  20050915 15:58:56

On Thu, 15 Sep 2005, Min Xu wrote: > It turns out the functionality requested _has_ been implemented in > MeshFunction. A while ago, Michael Schindler sent in a patch for MeshFunction > and other minor code corrections to this list. Has the patch been merged into > the official libmesh? I don't recall that particular patch, but it looks like MeshFunction is in the CVS libMesh at least. The mapping based point locator I described would be a little faster than the current PointLocator implementations, but it's an O(N) vs O(N log N) sort of difference. You should definitely start by using the existing code.  Roy 