From: KIRK, B. (JSC-E. (NASA) <ben...@na...> - 2005-02-04 17:49:19
|
Interesting application... Off the top of my head, I'd say solving the full Ke + Kb system simultaneously as one linear solve might be challenging. As you correctly point out, Kb is fully dense, while Ke would be a typical, sparse Laplacian discretization. The latter is amenable to preconditioned conjugate-gradient iterative techniques, while the former certainly would require a direct-type solver. Of course, you could just use a direct solver on the whole thing, but it would scale like O((N+n)^3), which is very bad for AMR. Enough rambling... From a library point-of-view I think we're in good shape. What I would suggest is *two* systems defined on two separate (but related) meshes. The first system would be your standard Laplace system on the volume mesh. The second system is the BEM system, defined on a 2D mesh in 3D space composed of the boundary faces of your volume mesh. You could then iterate a couple of times between the two systems until convergence. For example: 1.) Assume some initial boundary conditions, solve the Laplace system using these BCs. 2.) Solve the BEM system, using this Laplace system if needed in the forcing function. 3.) Re-solve the Laplace system with Neumann BCs given by the potential from the BEM solution. Repeat until convergence of each, separate system. It might sound wasteful, but I'm not sure it's that bad. For example, in step (1.) you build your sparse matrix Ke, form a preconditioner, and solve the linear system (with CG, probably). If you are using Neumann BCs & have a linear Laplace problem you can keep this system matrix & preconditioner, so you only form this small system once. Then, the bulk of the work in step (2.) will be forming the matrix and solving the dense system. Again, you might have to do this only once. For example, if you form an O(n^3) LU factorization in step (2.) then subsequent solves are only a forward and backward substitution sweep, which are O(n^2). In any case, this is certainly better than O((N+n)^3)! Finally, in step (3) of the iteration you probably only need to compute the RHS vector using the output of step (2). In particular, you may be able to reuse the matrix & preconditioner from (1) for this and subsequent solves. Not sure about what nonlinearities you have, but I would expect this approach to converge. It is like a non-overlapping Schwarz domain-decomposition scheme where you advance each subdomain in a Gauss-Siedel-like fashion. You could use mesh.create_submesh() and the boundary_side_iterators to build the 2D boundary mesh for the boundary system. This is probably what you need to do in either the approach I outlined or your proposed scheme. The only issue is that at the moment create_submesh() duplicates nodes rather than sharing them with the "parent" mesh. Some times you want this, but not in this case. I think I can get an easy fix together for that. In my scheme you would then create 2 systems, one for each mesh. For your scheme I'd need to smarten-up the DofMap to understand indexing multiple meshes. I've planned this anyway, a sort of CompositeMesh that contains several distinct meshes, but it's not yet implemented. Let me know what you think about this information overload... I think either approach might require a little change in the library, but not too much. I'd certainly be happy to help. -Ben -----Original Message----- From: lib...@li... [mailto:lib...@li...] On Behalf Of Mark Blome Sent: Friday, February 04, 2005 10:00 AM To: lib...@li... Subject: [Libmesh-users] hybrid FEM-BEM Hi everybody, I am currently implementing a way to deal with open boundaries in my FEM-solver for a laplace equation in 3D (well, lets say I am trying to do so). The idea is to use a hybrid FEM-BEM approach where the outer domain is calculated using a boundary element method (I use a variant of BEM called Trefftz method) while the inner domain is calculated using the conventional FEM galerkin formulation. The potential in the outer domain is approximated using a set of functions (which are fundamental solutions of the laplace equation), hence the DOFs in the outer domain are scalar factors to these functions. I would like to treat the outer domain as a "Superelement" with n DOFs associated to it. How can I handle the additional n DOFs ? How can I number these DOFs appropriately?The resulting equations should look something like: | Ke 0 | | u | = | Fe | | 0 Kb | | a | | Fb | where Ke, Fe are the conventional FEM stiffness matrix and RHS, Kb and Fb are the matrix and RHS for the BEM region (Kb is a nxn dense matrix, n=number of interpolating functions). To link the two regions together (the potential of the two solutions must equal on the boundary) I need to introduce constraint equations. They express each FEM DOF on the boundary by the n DOFs in "a" that represent the scalar factors to the functions I use for the interpolation of the potential field in the BEM region. I read about the constraint matrices that are used for AMR, but can I use these also for my case? If so, how do I use them? For example I would like to set these constraints locally for each element, but the DOFs in "a" are global. Thanks for any comments on this! Mark ------------------------------------------------------- This SF.Net email is sponsored by: IntelliVIEW -- Interactive Reporting Tool for open source databases. Create drag-&-drop reports. Save time by over 75%! Publish reports on the web. Export to DOC, XLS, RTF, etc. Download a FREE copy at http://www.intelliview.com/go/osdn_nl _______________________________________________ Libmesh-users mailing list Lib...@li... https://lists.sourceforge.net/lists/listinfo/libmesh-users |
From: Mark B. <mb...@au...> - 2005-02-07 09:57:41
Attachments:
hybrid_fem_trefftz_gyimesi.pdf
|
Hi Ben, this sounds like a very reasonable approach for a combined FEM-BEM. I read an article where they used Gauss-Seidel-Relaxation to assure the continuity of the potential but didn't really understand what they where talking about. In the approach that I would like to use, the exterior Region is approximated by a so called "complete T-Set" of functions (the method is called Trefftz). Rather than having one DOF per node like in conventional BEM there are only that much DOFs in the exterior region as T-Set functions are used to approximate the potential: U(r)=SUM( a_i * f_i (r) ), i=1..N The functions f_i (r) are chosen such that they a priori fullfill the laplace equation in the volume domain and hence only surface integrals have to be evaluated to estimate the unknowns "a_i". For accurate calculations something between 50 and 200 functions should be sufficient. Compared to the surface DOFs I usually have in my calculations (more than 5000) this is rather a small number and my hope is that this approach will allow me to solve the whole matrix in one step. So I would like to use the approach to solve one sytem in one step rather than solving two systems iteratively. If my idea will proof to be a bad one (in terms of speed) I will happily come back to your proposal! But as I am planning to use a direct solver for the whole system (cause I have to solve the same system for about 100 different RHS within each step of an data inversion code) I first want to try it that way (I've heard there are very nice direct solvers that make use of the sparsity pattern of a given matrix and hence should not be too wasteful with memory and cpu time). For my approach I need to extend the matrix with additional DOFs. I think the numbering should be quiet easy cause the Trefftz DOFs "a_i" may, as they are not linked to mesh nodes, just go to "the end of the matrix" and RHS. I guess the linking of the FEM boundary DOFs (potential values at nodes) to the DOFs "a_i" such that U_surfnode_fem = U_surfnode_trefftz <=> U_surfnode_fem = SUM( a_i * f_i (surfnode) ), i=1..N will be more difficult. Actually I have to say that I have no idea how to do this in libmesh. I only heard that there are constraint matrices. They seem to be appropriate for this, I guess ?!?! I will be happy to hear more comments from you on this, cheers, Mark PS: see the short article attached where the approach I am trying to use is described, just in case you want to read something about it.... |