Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

## Re: [Libmesh-users] libmesh and BEM

 Re: [Libmesh-users] libmesh and BEM From: Vasilis Vavourakis - 2011-01-11 16:34:52 ```2011/1/10 Roy Stogner > > On Mon, 10 Jan 2011, Vasilis Vavourakis wrote: > > (1) in BEM mostly boundary integrations occur, which means that in 2D >> problems only 1D discretization is needed. however, since only Edge >> elements >> are "loaded" to the MeshBase class, when i try to evaluate the normal >> vectors to each quadrature point for each boundary (edge) element through >> the "get_normals" function then i get an error. one option would be to >> load >> a 2D mesh (of triangles and quadrilaterals) and perform only boundary >> integration to all non-neighboring side elements, but this way i would >> load >> into the System the nodes located inside the domain. any ideas on that? >> > > The most flexible idea that jumps to my mind: adding a per-variable > flag to identify BEM variables in a System, then querying that in > DofMap so as to only assign degrees of freedom to boundary nodes. > This would be a difficult feature for a non-experienced libMesh > developer to write, though. > > i've thought of a simple but yet not smart enough solution: load two Meshes: one mesh that is 1D (this is the BEM mesh used for the system) and another dummy one that is 2D (with Trias+Quads+Edges), where there is one-to-one match on all Edge elements. in addition, force the same kind of order of Quadrature for both meshes for surface integrations...so that there is consistency to quadrature points. this way it can be evaluated the normal vectors...i think! BUT this solution hinders someone to use AMR directly to the 1D mesh without applying the same thing to the 2D mesh, if you understand what i mean. > > (2) the system matrix "A" of the final linear system "A*x=b" in not banded >> as in FEM, but full populated... this means that parallelization is not >> that >> easy to implement, right? >> > > Our SparseMatrix class would still parallelize large dense matrices; > you'd be able to get solutions. The only code change is again in > DofMap - the send_list for communication of a BEM variable would need > to be full rather than built up from neighbor topology. > right!!! there is no "neighbor topology" in BEM as in FEM ... the DofMap should address all nodes in the mesh!!! but about DoF index... should the function "dof_number()" do the job? for example: Node* node; ... const unsigned int index =node->dof_number(system.number(), system.variable_number("u"), 0); the above variable "index" is it local or global??? > Efficiency is another question, though; I'm not sure what you'd want > to do for preconditioning. > i think that direct solvers will do the job...so i've not come into thinking about pre-conditioners. the other problem that arises to this method is that the final system is * not* of the form: "K * u = F", as in FEM for elasticity; with "K" being the stiffness and "F" the external force vector. in BEM the final system is like: "H * u = G * t", where matrices "H", "G" contain the boundary integrals and "u", "t" is the displacement and traction vector to all boundary nodes. so in case there are some boundary nodes to whom the traction is unknown (known "u" : Dirichlet BC) then there is a need for some sort of system rearrangement!!! which to my mind seems very difficult for the moment to handle in parallel. Vas > --- > Roy > ```

### Thread view

 [Libmesh-users] libmesh and BEM From: Vasilis Vavourakis - 2011-01-10 17:50:57 ```hi all and happy new year :) i'd like to ask if anyone has implemented libmesh in Boundary Element Methods. i'm currently trying to built a simple 2D BEM solver for the Laplace PDE "Grad^2(u) = 0". since i'm not that a libmesh expert i face some difficulties in the implementation: (1) in BEM mostly boundary integrations occur, which means that in 2D problems only 1D discretization is needed. however, since only Edge elements are "loaded" to the MeshBase class, when i try to evaluate the normal vectors to each quadrature point for each boundary (edge) element through the "get_normals" function then i get an error. one option would be to load a 2D mesh (of triangles and quadrilaterals) and perform only boundary integration to all non-neighboring side elements, but this way i would load into the System the nodes located inside the domain. any ideas on that? (2) the system matrix "A" of the final linear system "A*x=b" in not banded as in FEM, but full populated... this means that parallelization is not that easy to implement, right? i have more questions but i think that the above are the most "difficult" ones i have so far... in case anyone in the past has worked in that area, their ideas and help would be much appreciated! many thanks, Vas ```
 Re: [Libmesh-users] libmesh and BEM From: Roy Stogner - 2011-01-10 18:52:46 ```On Mon, 10 Jan 2011, Vasilis Vavourakis wrote: > (1) in BEM mostly boundary integrations occur, which means that in 2D > problems only 1D discretization is needed. however, since only Edge elements > are "loaded" to the MeshBase class, when i try to evaluate the normal > vectors to each quadrature point for each boundary (edge) element through > the "get_normals" function then i get an error. one option would be to load > a 2D mesh (of triangles and quadrilaterals) and perform only boundary > integration to all non-neighboring side elements, but this way i would load > into the System the nodes located inside the domain. any ideas on that? The most flexible idea that jumps to my mind: adding a per-variable flag to identify BEM variables in a System, then querying that in DofMap so as to only assign degrees of freedom to boundary nodes. This would be a difficult feature for a non-experienced libMesh developer to write, though. > (2) the system matrix "A" of the final linear system "A*x=b" in not banded > as in FEM, but full populated... this means that parallelization is not that > easy to implement, right? Our SparseMatrix class would still parallelize large dense matrices; you'd be able to get solutions. The only code change is again in DofMap - the send_list for communication of a BEM variable would need to be full rather than built up from neighbor topology. Efficiency is another question, though; I'm not sure what you'd want to do for preconditioning. --- Roy ```
 Re: [Libmesh-users] libmesh and BEM From: Vasilis Vavourakis - 2011-01-11 16:34:52 ```2011/1/10 Roy Stogner > > On Mon, 10 Jan 2011, Vasilis Vavourakis wrote: > > (1) in BEM mostly boundary integrations occur, which means that in 2D >> problems only 1D discretization is needed. however, since only Edge >> elements >> are "loaded" to the MeshBase class, when i try to evaluate the normal >> vectors to each quadrature point for each boundary (edge) element through >> the "get_normals" function then i get an error. one option would be to >> load >> a 2D mesh (of triangles and quadrilaterals) and perform only boundary >> integration to all non-neighboring side elements, but this way i would >> load >> into the System the nodes located inside the domain. any ideas on that? >> > > The most flexible idea that jumps to my mind: adding a per-variable > flag to identify BEM variables in a System, then querying that in > DofMap so as to only assign degrees of freedom to boundary nodes. > This would be a difficult feature for a non-experienced libMesh > developer to write, though. > > i've thought of a simple but yet not smart enough solution: load two Meshes: one mesh that is 1D (this is the BEM mesh used for the system) and another dummy one that is 2D (with Trias+Quads+Edges), where there is one-to-one match on all Edge elements. in addition, force the same kind of order of Quadrature for both meshes for surface integrations...so that there is consistency to quadrature points. this way it can be evaluated the normal vectors...i think! BUT this solution hinders someone to use AMR directly to the 1D mesh without applying the same thing to the 2D mesh, if you understand what i mean. > > (2) the system matrix "A" of the final linear system "A*x=b" in not banded >> as in FEM, but full populated... this means that parallelization is not >> that >> easy to implement, right? >> > > Our SparseMatrix class would still parallelize large dense matrices; > you'd be able to get solutions. The only code change is again in > DofMap - the send_list for communication of a BEM variable would need > to be full rather than built up from neighbor topology. > right!!! there is no "neighbor topology" in BEM as in FEM ... the DofMap should address all nodes in the mesh!!! but about DoF index... should the function "dof_number()" do the job? for example: Node* node; ... const unsigned int index =node->dof_number(system.number(), system.variable_number("u"), 0); the above variable "index" is it local or global??? > Efficiency is another question, though; I'm not sure what you'd want > to do for preconditioning. > i think that direct solvers will do the job...so i've not come into thinking about pre-conditioners. the other problem that arises to this method is that the final system is * not* of the form: "K * u = F", as in FEM for elasticity; with "K" being the stiffness and "F" the external force vector. in BEM the final system is like: "H * u = G * t", where matrices "H", "G" contain the boundary integrals and "u", "t" is the displacement and traction vector to all boundary nodes. so in case there are some boundary nodes to whom the traction is unknown (known "u" : Dirichlet BC) then there is a need for some sort of system rearrangement!!! which to my mind seems very difficult for the moment to handle in parallel. Vas > --- > Roy > ```
 Re: [Libmesh-users] libmesh and BEM From: Roy Stogner - 2011-01-11 17:25:02 ```On Tue, 11 Jan 2011, Vasilis Vavourakis wrote: > load two Meshes: one mesh that is 1D (this is the BEM mesh used for the system) and another dummy one that is 2D (with Trias+Quads+Edges), where > there is one-to-one match on all Edge elements. in addition, force the same kind of order of Quadrature for both meshes for surface > integrations...so that there is consistency to quadrature points. this way it can be evaluated the normal vectors...i think! BUT this solution > hinders someone to use AMR directly to the 1D mesh without applying the same thing to the 2D mesh, if you understand what i mean.   Actually, if you're using a polygonal/polyhedral domain (or a domain where coarse piecewise quadratic geometry is good enough), the solution here is simple: don't refine the 2D mesh. You'll have to generate "quadrature points" via inverse_map rather than by hoping you matched up quadrature rules, but that's the safe way to go anyway (some of our 2D quadrature rules that you'd need on the boundary of a 3D domain aren't symmetric). A less hackish solution would be to just add (optional) tangent/normal vectors calculation for interior FE evaluations. The FEBase::compute_face_map code is straightforward enough to move and extend; we'd be happy to accept a patch. > Our SparseMatrix class would still parallelize large dense matrices; > you'd be able to get solutions.  The only code change is again in > DofMap - the send_list for communication of a BEM variable would need > to be full rather than built up from neighbor topology. > > > right!!! there is no "neighbor topology" in BEM as in FEM ... the DofMap should address all nodes in the mesh!!! > but about DoF index... should the function "dof_number()" do the job? for example: >   > Node* node; > ... > const unsigned int index =node->dof_number(system.number(), system.variable_number("u"), 0); > > the above variable "index" is it local or global??? The index is a global index. But getting indices isn't the problem; ensuring sufficient communication between every index is. That's what would require an additional System flag set and DofMap code. > the other problem that arises to this method is that the final system is not of the form: "K * u = F", as in FEM for elasticity; with "K" being > the stiffness and "F" the external force vector. in BEM the final system is like: "H * u = G * t", where matrices "H", "G" contain the boundary > integrals and "u", "t" is the displacement and traction vector to all boundary nodes. so in case there are some boundary nodes to whom the > traction is unknown (known "u" : Dirichlet BC) then there is a need for some sort of system rearrangement!!! which to my mind seems very difficult > for the moment to handle in parallel. This seems to assume that you're solving for u and t in the same function space, and that known coefficients for one correspond to unknowns for the other? In that case you just add_matrix and add_vector to allocate H,G,u_known,t_known terms, assemble them separately, copy the appropriate blocks into your final algebraic system to solve, then copy the results back after the solve. --- Roy Requires a little care, but there shouldn't be```