Periodic conditions in a weak form

[jfbarthelemy]

Dear Ferrite devs and users,

I saw that periodic boundary conditions can be defined on cell where nodes of a face are mirrors of nodes of its counterpart. But would it be possible to make it in a weak form (with Lagrange multipliers) in the case where the mesh is not defined with periodic information (not even the same number of nodes on periodic faces for instance)?

Thanks in advance

[fredrikekre]

There is no builtin support for this, but I suspect all the tools necessary to construct this on your own is available with some hacking. Depending on your formulation, I suspect you need an "integration mesh" and then map all the quadrature points there from reference coordinate -> spatial coordinate -> reference coordinate by looking up the correct element in the last step. This should be doable with the PointEvalHandler, I think.

[jfbarthelemy]

Thanks. Not sure that I would be able to do this now since I am too new to Ferrite but I keep in mind after some search in the tutorials and source.

[fredrikekre]

If you can come up with a minimal formulation (e.g. weak form) and post here I am sure someone will take the bite and post a solution.

[jfbarthelemy]

OK I consider a periodic unit cell such as that of the "computational_homogenization" tutorial except that the nodes of corresponding faces by periodicity (let us denote Γ⁺ and Γ⁻ a pair of corresponding faces) are not mirror from each other, I mean the mesh has been built without any mention of periodicity a priori but the boundary is of course still a square. The periodicity on the fluctuation displacement can be obtained by introducing a Lagrange mutliplier field λ (of vector type) on Γ⁻ and add the following term to the Lagrangian

$$ \int_{\Gamma^-} \lambda\cdot(u^+-u^-) dS $$

where $u^+$ is translated from Γ⁺ to the corresponding point of Γ⁻. Of course this should be done for the pair of right and left sides on the one hand and the pair of bottom and top sides on the other hand.

From a variational point of view this corresponds to the system

$$ \int_{\Omega} \delta\varepsilon:C:(\varepsilon^M+\varepsilon(u)) d\Omega= \int_{\Gamma^-} \lambda\cdot(\delta u^+-\delta u^-) dS $$

$$ \int_{\Gamma^-} \delta\lambda\cdot(u^+-u^-) dS=0 $$

where the unknown fields are u (defined on the whole domain) and $\lambda$ (defined only on Γ⁻).

To be able to handle such a problem, one should be able to

1. build a variable living only on a part of the boundary (ie faceset and not cellset)
2. translate the field u from one face to the other face to perform the integration on only one of the mirror faces

This is typically what is done in GetFEM with interpolate transformations here.

I don't know if this is clear enough. Don't hesitate to tell me if it isn't so.

Let me also precise that although I am aware that periodic conditions can be better handled by the method provided in "computational_homogenization" tutorial, I also have other examples where I would find relevant to have a variable field living only on a faceset on the boundary but also possibly inside the domain. As an example, I would mention the problem of electrochemistry: I have on ohmic material in the bulk domain (typically concrete) and a metal layer on a part of the boundary. The metal is much more conductive than concrete but the fact that its geometry corresponds to a thin layer, it is modeled as a faceset with a surface conductivity. The variables are then the electric potential within the concrete domain and electric potential on the metal facesets and there exists a relationship between the current density from metal to concrete and the local potential jump. This problem is well posed (with other classical boundary conditions of Dirichlet and/or Neumann types elsewhere of course) but it requires mixed multivariables. Would it be possible with Ferrite, now or in the future?

Thanks

[termi-official]

I think this is possible with the following approach using the current master, which might require some work. We do not have direct support for trace spaces in Ferrite.jl directly.

1. Let us say you have a 3D problem with hexahedral elements, then you can add helper quadrilateral elements on the boundary of interest and use the SubDofHandler to define the same fields on the Hex element subdomain+Quad element "boundary subdomain". This should give you enough control for this part. Just note that for the next step it is necessary to add the "face quads" in a way such that they match the corresponding orientation of the Hex thy are attached to.
2. I am not sure if there is a way to automate the translation, including finding the translation function for a general geometry, but assuming you can define a translation function from one boundary to another one, then you should be able to get away by translating the integration points on one boundary in physical coordinates (not reference). Then the "only" missing part is the transformation from physical to reference on the opposing element. For linear geometries there is the The evaluation can be done vie PointEvalHandler (this will fail for curved geometries, because it is just not implemented, yet). However, it currently does not work for the face elements directly. So, what you can do is to find the quadrature point of the Hex to which the face is attached to with the PointEvalHandler. Now with the quadrature point given at the Hex you can transform it by building the inverse of this function here https://github.com/Ferrite-FEM/Ferrite.jl/blob/master/src/FEValues/face_integrals.jl#L128-L137 (which currently transforms from a given face to an element).

Does this give you an idea on how this can be implemented?

[jfbarthelemy]

Thanks. I think it is clear but I am not sure that I will get into such a tedious workflow now as I am discovering Ferrite. I will stick to what is already simply available and stay informed for future improvements. This discussion made me better understand the code which is already very nice.