Using the solution from one problem as a boundary condition in another problem with matching mesh on the boundary

[krei]

The problem I am currently having is directly related to my previous post [here](https://groups.google.com/forum/#!searchin/dealii/krei/dealii/eq_zP0jrSJU/SyQXA7A-DAAJ).

The physics are exactly the same as described in that thread: I have a two domain system: vacuum and metal, both domains have their own meshes. [Here is a picture](http://i.imgur.com/CejXi1y.png). The meshes are exactly matching on the boundary. I want to solve the Laplace equation for electric fields in the vacuum part (does not depend on anything in the metal part) and then use the electric field as a boundary condition to the metal part.

I am currently using VectorTools::point_gradient to evaluate the electric field from the Laplace problem in the quadrature points of the boundary cell faces of the metal mesh. As a result, it takes over an hour to assemble the system (solving takes less than a second). 

Now, considering that the mesh is exactly matching at the boundary (all the quadrature points etc), how could I efficiently evaluate the electric field at the boundary?

And I apologize if I'm mistaken, but there doesn't seem to be any tutorials on these kinds of problems where you have multiple domains with multiple meshes. Or is it usually done by using a single mesh and somehow defining different domains in that mesh? (haven't delved deeply into Step-46, but there they do something like this, right?)

[Daniel Arndt]

krei,

If you want to solve different PDEs on different domains that can be discretized by a common mesh, the preferred approach is to use a hp-vector finite element.

This means that on each of your subdomains all blocks of your finite element but one are of type FENothing.

You might want to have a look into step-46 for how to do this.

Best,

Daniel

[krei]

Thanks for the response. I think this is worth a try. However, as the electric field calculation is entirely decoupled from the nonlinear currents & heating part, which need to be solved by newton's method, then wouldn't it be effective to just once calculate the field and then start the newton iterations for only temperature and currents (probably the performance can be made the same, but perhaps code readability decreases)? And the second thing is that I want a modular code, where you can choose if you want to only calculate electric fields and disregard currents & heating or calculate both, in this case it would also be good to keep the codes separate. So it would be really convenient to have a some kind a fast function like "evaluate_gradient_in_quadrature_point(int index)" or something similar.

[krei]

Hello,

I mostly implemented the hp-vector finite element approach (according to step-46), but alas, I think it might not be applicable. (I simplified the boundary condition in original post a bit.) In my case I need to apply an emission current boundary condition to the electric currents in copper based on the electric field in vacuum. The emission currents are not expressed through the electric field by a simple manner, more specifically, J ~ E^2, J ~ 1/E and J ~ exp(a*E)  and I might want to use interpolation from a grid to evaluate J(E). (J - emission current density, E - electric field at the surface).

If I want to solve everything in one big system (as is done in step-46), then I don't have access to the electric field values during the system assembly, I can access shape functions but I can't evaluate the emission current through them.

Perhaps I have misunderstood and I can somehow evaluate the boundary conditions? Or what would a better approach be? 

On Saturday, July 30, 2016 at 1:38:06 PM UTC+3, Daniel Arndt wrote:

[Daniel Arndt]

krei,

If your emission current boundary conditions do not depend linearly on the electric field, the whole problem becomes non-linear and hence you can't solve the whole problem directly.
What you can do is to first solve for the electric field and afterwards for the metal part. In particular, this means to neglect in the first assembly all the cells that a related to the metal part.
After solving you would then loop over the interface to find out which dofs on the metal part you want to constrain. For this you would initialize a Quadrature object with the unit support_points [1].
At the same time you can use FEValues::get_function_* to evaluate the electric field on the locations of these dofs. This should be much faster than VectorTools::PointGradient.
Finally, assemble and solve for the metal part taking the computed constraints into account.

Does this make sense to you?

Best,
Daniel

[1] https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element

[Michael Harmon]

Hey,

I had a similar problem: PDES in separate domains that are coupled through an interface as a boundary condition.  You can go about it using one triangulation; I attempted to do this at first, but ended up using multiple meshes. The fact you have matching meshes on the boundary is good.  What you want to do is create two meshes and two dof handlers.  Initially you need to create a std::map to create a bijection between the faces on the different meshes that are on the interface/ connecting boundary.  When you are assembling the coupling terms on the boundary condition you can get the corresponding face in the other mesh and initialize the fe_values on that face and get its values on the face to assemble the boundary conditions.

Check out:

https://groups.google.com/forum/#!searchin/dealii/multiple$20triangulations%7Csort:relevance/dealii/R64BUBhyFf0/LbbxH9v8GgAJ

If it helps I can also send you my code to see how to do this using deal.ii

Best,

Mike

[krei]

Thanks for the response. I have a more-or-less working version where the whole thing is solved together by Newton's method, but I bet solving electric field separately gives a good performance boost. I will try to fiddle with it, however, generating multiple meshes with matching boundary for testing is not very trivial (GMSH doesn't let you do this, as far as I know), I have to investigate if I can mark the elements in gmsh and then separate the meshes in deal ii.

[Michael Harmon]

Nice! Good luck with Gmsh!

[vladisla...@gmail.com]

Hello, Mike,

I have a similar problem like you. I have two PDEs in separate domains coupled through an interface. I have also tried to solve this using step-46.However, my code haven't produced the correct solution in my test case. Thus, I am very interested to try new approaches, especially how you have solved the problem using two different meshes (which of course fit on the interface). Could you send me a code snippet showing how you deal with the interface term and the vertices on the different meshes? It would be very helpful. Your idea using std::map makes sense to me, however I still can not really imagine how everything works together in the code.

Best, 
Vladislav

[Kristjan Eimre]

Hi Vladislav,

You can take a look at my project here: https://github.com/eimrek/dealii-field-currents-heating, perhaps it's useful. In currents_and_heating.h/.cc files I use a map between vacuum and metal face cells to evaluate boundary condition from one domain in the other. Note that in my case the coupling is only one-directional.

Mike's repository is here: https://github.com/mdh266/PECS, this might be also useful for you.

Best,
Kristjan

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/Uvi3-xhfHBk/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[vladisla...@gmail.com]

Thank you, Kristjan.

Both links are helpful.  

To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.