FSI driver

[fverdugo]

Hi @santiagobadia @oriolcg,

In order to test the new surface-coupling machinery, the best way is to write a simple FSI tutorial. To this end I have created a new empty tutorial here https://github.com/gridap/Tutorials/blob/fsi_tutorial/src/fsi_tutorial.jl

It is in branch fsi_tutorial.

Could you design a simple FSI problem, and chose a simple FSI formulation? You can start writing the equations in the file above if you like.

Once the problem and formulation are specified I can tell if it can be implemented with the current machinery, and then we can start the implementation.

Do you want to restrict to a steady-state problem for simplicity or do you prefer a transient one?

[santiagobadia]

There are two types of numerical surface coupling we want to consider.

1. Globally defined unknowns

The first type of coupling is the one in which we have different problems in different regions, continuity is enforced via unknowns that are enforced to belong to a global FE space in the whole domain. Continuity of traces is straightforwardly imposed via the continuity of the FE space on the whole domain. Continuity of fluxes is naturally imposed after integration by parts of the different problems at hand, eliminating the interface Neumann-type terms. I would start with this one. In order to check that everything works, I would use the vector Laplacian coupled to the Stokes problem. It is the simplest problem I can imagine that is related to FSI. To change the vector Laplacian with linear elasticity in irreducible form is a question of minutes. I would not consider transient problems or moving meshes yet. It can be done once the approach is proven successful.

2. Weak coupling

The second approach is to define different FE spaces at the different regions, and enforce the continuity in weak form using Nitsche or DG-like terms on interfaces. We can use the same test as above, consider two different spaces for the velocity (one in each region) and enforce continuity on the interface adding these interface terms.

Other approaches that we can consider in the future...

3. Strong coupling through constraints

I guess there is room for other options, e.g., to enforce the continuity strongly but when constraints are not us = uf but something more complicated, e.g., d1 = dt*uf+d0. It can probably be done through Dirichlet boundary conditions... to think about it.

4. Hybrid formulations

Another approach is to define a space of traces on a mesh for the interface, univalued, and enforce continuity on the other two sides in terms of this one...

[oriolcg]

I can take care of building this tutorial, following the steps that @santiagobadia outlined in the previous comment.

[santiagobadia]

I see that 1. is already in this test.

Thus, for 1., the answer is Yes, and it is working. A very minor thing is the need to introduce mean operator for a function v that is continuous... Probably a face_restrict to make it clearer. However, I don't think that would be needed, see e.g. https://github.com/gridap/Tutorials/blob/5fd09b689aba7e153edb3e887d0f6be5520de8f2/src/t005_dg_discretization.jl#L184

@oriolcg you could probably take this tutorial as starting point, replace the global velocity space by two indepedent spaces in the two different solid and fluid regions, and enforce continuity weakly using the terms as in the DG tutorials. What do you think? My only concern is the jump operator that will be needed there. I don't think it will work with the current machinery. I would write the terms explicitly without jump, BUT I am not sure the mean operator could be used there (functions only defined in one side, not sure it will return an error or 0.5 the value). If mean does not work (even if it works, it would be misleading), a face_restrict operator (or nothing) would be all we need there. I think it is very simple, just a simplification of mean.

But let us wait what @fverdugo has in mind first.

[santiagobadia]

I have filled the tutorial with a compressible linear elasticity + Stokes (incompressible linear elasticity) problem. @oriolcg you could start from here, to fill the tutorial, etc.

[oriolcg]

I have updated the problem formulation, but it looks like I don't have permission to push changes to this repository. How should I proceed?

I'll leave the formulation here for you to have it in the meantime

# ## Problem statement
#  Let $\Gamma_{\rm FS}$ be the interface between a fluid domain $\Omega_{\rm F}$ and a solid domain $\Omega_{\rm S}$. We denote by $\Gamma_{\rm F,D}$ and $\Gamma_{\rm F,N}$ the Dirichlet and Neumann conditions on the fluid boundary.
# The Fluid-Structure Interaction (FSI) problem reads: find $ u_{\rm F} $, $ p_{\rm F} $ and $ u_{\rm S} such that
# ```math
# \left\lbrace
# \begin{aligned}
# -\nabla\cdot\boldsymbol{\sigma}_{\rm F} = f &\text{ in }\Omega_{\rm F},\\
# \nabla\cdot u_{\rm F} = 0 &\text{ in } \Omega_{\rm F},\\
# -\nabla\cdot\boldsymbol{\sigma}_{\rm S} = s &\text{ in }\Omega_{\rm S},\\
# \end{aligned}
# \right.
# ```
#
# satisfying the Dirichlet and Neumann boundary conditions
# ```math
# \left\lbrace
# \begin{aligned}
# # u_{\rm F} = g &\text{ on } \Gamma_{\rm F,D},\\
# \boldsymbol{\sigma}_{\rm F}\cdot n_{\rm F} = 0 &\text{ on } \Gamma_{\rm F,N},\\
# \end{aligned}
# \right.
#```
#
# and the kinematic and dynamic conditions at the fluid-solid interface
# ```math
#  u_{\rm F} = u_{\rm S} &\text{ on } \Gamma_{\rm FS},\\
# \boldsymbol{\sigma}_{\rm F}\cdot n_{\rm F} + \boldsymbol{\sigma}_{\rm S}\cdot n_{\rm S} = 0 &\text{ on } \Gamma_{\rm FS}.\\
# \end{aligned}
# \right.
#```
#
# Where $\boldsymbol{\sigma}_{\rm F}=2\mu_{\rm F}\boldsymbol{\varepsilon}(u_{\rm F}) - p_{\rm F}\mathbf{I}$ and $\boldsymbol{\sigma}_{\rm F}=2\mu_{\rm S}\boldsymbol{\varepsilon}(u_{\rm S}) +\lambda_{\rm S}tr(\boldsymbol{\varepsilon}(u_{\rm S}))\mathbf{I}$.

[santiagobadia]

@fverdugo please find a Nitsche tutorial in which I have put a couple of comments about unclear things.

Problem 1, restrictions of FE spaces with strong Dirichlet data, triangulation vs model, etc.

https://github.com/gridap/Tutorials/blob/21ae7d91ed2974641de3b974698f14625db4aa76/src/fsi_tutorial_nitsche.jl#L97-L106

Problem 2,

Can we integrate the Nitsche coupling terms https://github.com/gridap/Tutorials/blob/b093ed83004c8cd1310f51dcdaf2081838a55340/src/fsi_tutorial_nitsche.jl#L189-L196 using https://github.com/gridap/Tutorials/blob/b093ed83004c8cd1310f51dcdaf2081838a55340/src/fsi_tutorial_nitsche.jl#L216-L217 ?

@oriolcg now you have write permits

[fverdugo]

What I have implemented so-far is the strategy "1. Globally defined unknowns" that @santiagobadia mentioned. I am also aware of the other strategies, but for now I have considered this as a first step.

Imposing interface continuity a la nitsche is the tricky part still not supported elegantly, thus I would not recommend to start writing a tutorial that uses this formulation. The key is that the interface operators jump and mean only apply to a single field.... A possible option is to define a new global space for fluid and solid that is continuous within the sub-regions and discontinuous at the interface

Something like:

 V = TestFESpace(  
   model=model, 
   valuetype=VectorValue{2,Float64}, 
   reffe=:QLagrangian, 
   order=order, 
   conformity =:H1, 
   discontinuous_at = trian_Γ, # Note this argument (not yet implemented)
   dirichlet_tags="dirichlet") 

Then, we can use jump and mean as always. By the way, it is also possible to apply a different weight to each side of the interface if necessary, thus I think this is quite general.

The second option (in fact the two options can co-exist) is to define two different (extended) spaces for fluid and another for solid. Currently it is possible to view the interface triangulation either from the left or the right side:

trian_Γ = InterfaceTriangulation(model,...)
trian_Γ_fluid = get_left_boundary(trian_Γ)
trian_Γ_solid = get_right_boundary(trian_Γ)

The result trian_Γ_fluid and trian_Γ_solid are standard boundary triangulations like the ones we already use for neumann BCs

If uh_fluid and uh_solid are quantities that are !=0 in the corresponding region and ==0 outside (i.e., extended quantities) one can always do:

uh_Γ_fluid = restrict(uh_fluid,trian_Γ_fluid)
uh_Γ_solid = restrict(uh_solid,trian_Γ_solid)

and then implement the terms in the weak form in function of uh_Γ_fluid and uh_Γ_solid as needed.

The only minor missing part is a new type of FETerm that allows the user to restrict each field (solid/fluid) to a different triangulation (now all fields are restricted internally to the same triangulation). This is not a big deal.

Strategy 3: to-think.

[santiagobadia]

@fverdugo my point is that you do not need mean and jump for surface coupling, see how I would write the Nitsche coupling. I think that this is more than acceptable. I think that creating multi-field operators for different unknowns is not even how people write it in papers (I should check this last statement, but I think it is generally true). In fact, you will want to consider different terms at different sides, related to jumps of coefficients, the fact that you are coupling a velocity and a displacement, etc.

[fverdugo]

Yes, this is the second option I commented above.

However, this does not work:

trian_Γ = InterfaceTriangulation(model,...)
quad_Γ = CellQuadrature(trian_Γ,degree)
 function nitsche_Γ(x,y) 
   us, uf, p = x 
   vs, vf, q = y 
   (γ/h)*vf*uf - vf*(n_Γ*ε(uf)) - (n_Γ*ε(vf))*uf + (p*n_Γ)*vf + (q*n_Γ)*uf 
      - (γ/h)*vf*us + (n_Γ*ε(vf))*us - (q*n_Γ)*us 
      + (γ/h)*vs*us + vs*(n_Γ*ε(us)) + (n_Γ*ε(vs))*us 
      - (γ/h)*vf*uf - (n_Γ*ε(vf))*uf + (q*n_Γ)*uf 
 end 
t_Γn_nitsche = FESource(nitsche_Γ,trian_Γ,quad_Γ)

~The good new is that I have just realized that this should work:~

trian_Γ = InterfaceTriangulation(model,...)
quad_Γ = CellQuadrature(trian_Γ,degree)
 function nitsche_Γ(x,y)

  # Take left or right side depending on the case
   _, uf, p = x.left
  us, _, _ = x.right 
   _, vf, q = y.left
   vs, _, _ = y.right

   (γ/h)*vf*uf - vf*(n_Γ*ε(uf)) - (n_Γ*ε(vf))*uf + (p*n_Γ)*vf + (q*n_Γ)*uf 
      - (γ/h)*vf*us + (n_Γ*ε(vf))*us - (q*n_Γ)*us 
      + (γ/h)*vs*us + vs*(n_Γ*ε(us)) + (n_Γ*ε(vs))*us 
      - (γ/h)*vf*uf - (n_Γ*ε(vf))*uf + (q*n_Γ)*uf 
 end 
t_Γn_nitsche = FESource(nitsche_Γ,trian_Γ,quad_Γ)

~I think that this API is quite nice and it should already work. Is this what you need?~

[fverdugo]

Sorry is the other way around (first unpack the fields and then take the left/right):

 function nitsche_Γ(x_Γ,y_Γ)

  us_Γ, uf_Γ, p_Γ = x_Γ
  uf = uf_Γ.left
  p = p_Γ.left
  us = us_Γ.right

  vs_Γ, vf_Γ, q_Γ = y_Γ
  vf = vf_Γ.left
  q = q_Γ.left
  vs = vs_Γ.right

   (γ/h)*vf*uf - vf*(n_Γ*ε(uf)) - (n_Γ*ε(vf))*uf + (p*n_Γ)*vf + (q*n_Γ)*uf 
      - (γ/h)*vf*us + (n_Γ*ε(vf))*us - (q*n_Γ)*us 
      + (γ/h)*vs*us + vs*(n_Γ*ε(us)) + (n_Γ*ε(vs))*us 
      - (γ/h)*vf*uf - (n_Γ*ε(vf))*uf + (q*n_Γ)*uf 
 end 

[fverdugo]

The following API is not implemented but it would take me 15 mins literally.

function nitsche_Γ(x,y)

  # Take left or right side depending on the case
   _, uf, p = get_left_side(x)
  us, _, _ = get_right_side(x) 
   _, vf, q = get_left_side(y)
   vs, _, _ = get_right_side(y)

   (γ/h)*vf*uf - vf*(n_Γ*ε(uf)) - (n_Γ*ε(vf))*uf + (p*n_Γ)*vf + (q*n_Γ)*uf 
      - (γ/h)*vf*us + (n_Γ*ε(vf))*us - (q*n_Γ)*us 
      + (γ/h)*vs*us + vs*(n_Γ*ε(us)) + (n_Γ*ε(vs))*us 
      - (γ/h)*vf*uf - (n_Γ*ε(vf))*uf + (q*n_Γ)*uf 
 end 

[santiagobadia]

I think I understand the problem. Even though us only lives in one side of Γ (in the other side it is conceptually zero, but it has no associated DoFs, etc), we need to extract the correct non-zero side. This is not needed on the domain boundary, because it is uniquely defined. Am I right?

However, I am not sure I understand why us requires right and uf requires left. It could be the other way around, couldn't it? How the user can ascertain which is the right case?

[fverdugo]

I think I understand the problem. Even though us only lives in one side of Γ (in the other side it is conceptually zero, but it has no associated DoFs, etc), we need to extract the correct non-zero side. This is not needed on the domain boundary, because it is uniquely defined. Am I right?

That's it

However, I am not sure I understand why us requires right and uf requires left. It could be the other way around, couldn't it? How the user can ascertain which is the right case?

Yes of course, it depends on how the InterfaceTriangulation(model, cell_to_is_fluid) was constructed. By convention, cells with a true value in the boolean array cell_to_is_fluid are on the left side. Thus, fluid is on the left in this example. Just a convention, but one needs some kind of convention.

[fverdugo]

Another variant that should work right now:

 function nitsche_Γ(x,y)

  uf = x[2].left
  p = x[3].left
  us = x[1].right

  vf = y[2].left
  q = y[3].left
  vs = y[1].right

   (γ/h)*vf*uf - vf*(n_Γ*ε(uf)) - (n_Γ*ε(vf))*uf + (p*n_Γ)*vf + (q*n_Γ)*uf 
      - (γ/h)*vf*us + (n_Γ*ε(vf))*us - (q*n_Γ)*us 
      + (γ/h)*vs*us + vs*(n_Γ*ε(us)) + (n_Γ*ε(vs))*us 
      - (γ/h)*vf*uf - (n_Γ*ε(vf))*uf + (q*n_Γ)*uf 
 end 

[fverdugo]

@santiagobadia

which API do you prefere?

For me, this is a good one and its already working:

function nitsche_Γ(x,y)

  uf = x[2].left
  p = x[3].left
  us = x[1].right

  vf = y[2].left
  q = y[3].left
  vs = y[1].right

   (γ/h)*vf*uf - vf*(n_Γ*ε(uf)) - (n_Γ*ε(vf))*uf + (p*n_Γ)*vf + (q*n_Γ)*uf 
      - (γ/h)*vf*us + (n_Γ*ε(vf))*us - (q*n_Γ)*us 
      + (γ/h)*vs*us + vs*(n_Γ*ε(us)) + (n_Γ*ε(vs))*us 
      - (γ/h)*vf*uf - (n_Γ*ε(vf))*uf + (q*n_Γ)*uf 
 end 

[fverdugo]

Another convention is that

trian_Γ = InterfaceTriangulation(model, cell_to_is_fluid)
n_Γ = get_normal_vector(trian_Γ)

returns the normal outwards to the left side.

[santiagobadia]

Ok this one

 us_Γ, uf_Γ, p_Γ = x_Γ
  uf = uf_Γ.left
  p = p_Γ.left
  us = us_Γ.right

or the one above (it is the same mechanism). Not sure we could think about something more expressive than right and left, but I agree that we need a convention, inward and outward, I don't know... minor thing.

[santiagobadia]

By the way, take a look at the other problem, the fact that we need a model restriction for FE spaces with strong Dirichlet data. If we can solve the problem we have Nitsche working. It is related to the other issue I have opened in Gridap.

[fverdugo]

By the way, take a look at the other problem, the fact that we need a model restriction for FE spaces with strong Dirichlet data. If we can solve the problem we have Nitsche working. It is related to the other issue I have opened in Gridap

Already implemeted (perhaps API can be further improved):

model = CartesianDiscreteModel(...)

model_solid = DiscreteModelPortion(model, list_of_solid_cells)
model_fluid = DiscreteModelPortion(model, list_of_fluid_cells)
...
trian = Triangulation(model)
trian_solid = RestrictedTriangulation(trian, cell_to_is_solid)
trian_fluid = RestrictedTriangulation(trian, cell_to_is_fluid)

V_f = TestFESpace(
  model=model_fluid,
  valuetype=VectorValue{D,Float64},
  order=order,
  reffe=:Lagrangian,
  conformity=:H1,
 dirichlet_tags = [1,2,5],
  restricted_at=trian_fluid)

V_s  = TestFESpace(
  model=model_solid,
  valuetype=VectorValue{D,Float64},
  order=order,
  reffe=:Lagrangian,
  conformity=:H1,
 dirichlet_tags = "boundary",
  restricted_at=trian_solid)

The restricted_at argument signals that the model is defined in a sub-region and that you want to extend the resulting FE space to the whole domain with a 0-extension (what you want for surface-coupling). If you don't pass restricted_at, the resulting FESpace will be defined only in the sub-region (which is what we need for embedded FEM)

A possible way to improve the API is to remove the restricted_at and define a new type RestrictedDiscreteModel that would signal that you want to perform the 0-extension.

We could end up having:

• DiscreteModelPortion / TriangulationPortion (both already implemented) for the cases in which you don't want to perform the 0-extension (i.e., for embedded FEM)
• RestricedDiscreteModel (not implemented) / RestricedTriangulation (already implemented) for the cases in which you want to perform the 0-extension (i.e. surface coupled multi-physics)

then the restricted_at argument will not be needed anymore

[fverdugo]

@oriolcg @santiagobadia today a lot of changes have been introduced in Gridap#master with respect to surface-coupling. See the issue https://github.com/gridap/Gridap.jl/issues/204 and the test https://github.com/gridap/Gridap.jl/blob/master/test/GridapTests/SurfaceCouplingTests.jl for more info

[fverdugo]

I have taken a look to the two fsi tutorials we have in this branch. I see that you are considering a toy problem with manufactured solution. For the tutorials, I usually prefer to consider problems with a little bit more of engineering meaning (test with manufactured solutions are already in Gridap/test). In this case, we could consider for instance a flow around a flexible flag:

For the moment steady state, and small displacements for the solid. Even though it would be also nice to use an ALE formulation. @santiagobadia do you believe we can already implement an ALE with the current machinery?

I also see that we are using the CartesianDiscreteModel. However, it would be also nice to use a model generated with GMSH in which we define the "solid" and "fluid" tags since this is what is done in practice.

[oriolcg]

I updated the fsi_tutorial with the problem that you propose. Here I used the cartesian mesh as a first step, I'll create a GMSH model in following updates.

Regarding the ALE formulation, I don't know if it makes much sense in a steady problem. Don't you think we should have an unsteady version of the tutorial first?

[santiagobadia]

ALE for steady is not possible. Another story is to consider a nonlinear steady problem in which the fluid domain changes due to the solid domain (assuming large displacements). But I think it does not pay the price to do this. Better wait till transient problem, I have already talked about that with @fverdugo .

[santiagobadia]

What can be done next is to write a problem coupled via Nitsche, using the machinery @fverdugo has already implemented. As far I can see, we are ready to implement the coupled FSI Nitsche problem as discussed in one of the posts above. I am missing something @fverdugo ?

[fverdugo]

Nice pictures @oriolcg. Next step is to solve a more complex geometry with GMSH, e.g., flow around a COVID-19 :-)

I really didn't mean ALE since it is a steady problem. I was thinking in defining a displacement in the fluid domain (without any feedback on the fluid) just as a first step for having an ALE in unsteady cases. It can be still a linear problem if we assume small displacements. In this scenario, it is perhaps better to define solid displacement + "fluid displacement" in the same field everywhere, and fluid velocity and pressure in two other fields restricted to the fluid domain. In any case, I would not include this for the moment.

For Nitsche, when integrating on the interface you only need to take the values at the corresponing side of the interface (left/right or inward/outward) as we have commented in some post above. It should work, even for non-oriented meshes.

[fverdugo]

BTW @oriolcg, with the latest version in Gridap#master, you can replace RestrictedTriangulation with Triangulation and remove the restricted_at key-word argument in the TestFESpace constructor

[oriolcg]

I started adding the approach with different FE spaces for fluid and solid velocities (see new commit). At the moment I'm merging both approaches in the same tutorial (case A: a unique FE space, case B: a fluid and a solid velocity FE spaces). Once everything works we can keep both approaches, split them in different tutorials or remove case A.

In principle I adapted everything, but I'm getting an assertion at line https://github.com/gridap/Gridap.jl/blob/08a30f171cf62ccbbda808fe83282531397958fd/src/FESpaces/FESpaceFactories.jl#L196 I think that the option of using the triangulation instead of a model is only implemented for L2 conformity, not for H1. Am I wrong?

I also tried to replace RestrictedTriangulation by Triangulation and I get an error. I'll try to go deeper and see if I can figure out what's happening.

[fverdugo]

I also tried to replace RestrictedTriangulation by Triangulation and I get an error. I'll try to go deeper and see if I can figure out what's happening.

Make sure you are using the latest Gridap version in master

(Tutorials) pkg> add Gridap#master

I think that the option of using the triangulation instead of a model is only implemented for L2 conformity, not for H1. Am I wrong?

Yes, continuous spaces require all the topological information provided by models, a triangulation (i.e. a simple integration mesh) is not enough.

By the way, with the newest version in Gridap#master you can easily create the model for the fluid part:

# choose the one fits better to your case
model_fluid = DiscreteModel(model,cell_to_isfluid)
model_fluid = DiscreteModel(model,list_of_fluid_cells)
model_fluid = DiscreteModel(model [, labels], "fluid")

moreover, if you have a model for the solid and another for the fluid, you can build the interface simply as

trian_fsi = InterfaceTriangulation(model_fluid,model_solid) 

[oriolcg]

Hi @fverdugo , I pushed a couple of commits in fsi_tutorial branch with further developments. However, I still have a couple of errors that I'm struggling to debug. I created a test that replicates the errors that I get in the tutorial, see https://github.com/gridap/Tutorials/blob/fsi_tutorial/test/fsi_failing_test.jl

The first error that I get is:

ERROR: LoadError: MethodError: no method matching symmetic_part(::Float64)

Uncommenting line 112 and commenting line 111 on the fsi_failing_test.jl I get the following error:

ERROR: LoadError: MethodError: no method matching array_cache(::Gridap.Geometry.SkeletonPair{Gridap.MultiField. ...

Any suggestions on how to solve them?

[santiagobadia]

The first one seems easy, symmetic -> symmetric :-)

[santiagobadia]

I am not sure I understand the second problem. If you comment 111 and uncomment 112, you are only integrating a right-hand side, no matrix. So, I expect it to fail when trying to solve the system...

[oriolcg]

Good catch with the symmetic typo @santiagobadia. But this is deep inside Gridap, in src/TensorValues/Operations.jl and src/Fields/FieldOperations.jl. I think the error is related to the fact that it is receiving a float instead of a tensor, but I cannot identify the origin of this error.

Regarding the second error, I get the same in the tutorial where we do have a LHS even discarding the Nitsche terms. Do you think the error is still caused by the fact that I'm commenting the Nitsche terms in the tutorial?

[santiagobadia]

OK, it seems that symmetic is in Gridap. I will fix it. But the error is simple, it cannot find the method symmetic_part. In which subroutine does it appear?

With regard to the second one, https://github.com/gridap/Tutorials/blob/0d2dac3089c9a32f86f20c2b3e80f931b190ac7f/test/fsi_failing_test.jl#L106-L109

you define a RHS only. So, if you do https://github.com/gridap/Tutorials/blob/0d2dac3089c9a32f86f20c2b3e80f931b190ac7f/test/fsi_failing_test.jl#L112-L113 it cannot work, you have no LHS. Where are the rest of the terms? By the way, where are the interior cell terms?

[oriolcg]

I created test/fsi_failing_test.jl with the minimum code required to replicate the error that I was getting in src/fsi_tutorial.jl, just to simplify the debugging process. All the terms are in the tutorial: https://github.com/gridap/Tutorials/blob/0d2dac3089c9a32f86f20c2b3e80f931b190ac7f/src/fsi_tutorial.jl#L261-L317

[fverdugo]

OK, it seems that symmetic is in Gridap. I will fix it. But the error is simple, it cannot find the method symmetic_part. In which subroutine does it appear?

Sorry for the typos... I am very prone to do typos, and the spell checker does not work for variable names...

When fixing these kind of typos, if the fixed name is exported in Gridap or in some submodule, I would recommend to define a constant with the old name:

const symmetic_part = symmetric_part

This allows us to publish a backwards compatible version of the library with the fixed typos, e.g., now the new version would be 0.8.1. Then, in the next backwards incompatible version, i.e., 0.9.0, we can remove these constants.

[fverdugo]

@oriolcg @santiagobadia I have just realized that the current FETerm machinery does not support to take the left/right side of the boundary as we are trying to do. In theory, writing jump(u_Γ) or -jump(u_Γ) instead of u_Γ.left should work unless we have a bug (if I am not wrong jump = left - right) . For gradients, symmetric gradients, etc: jump(gradient(u_Γ)) instead of gradient(u_Γ.left).

I have tried this in your failing test, and still does not work but now it seems to be a much minor bug.

Of course, we also want a way of writing these kind of terms without using the jump function. I am thinking perhaps in creating an API like:

function a(x,y)
# these would be  already in the left / rigth side
us,uf,p = x
vs,vf,q = y
# your terms here
end

t = LinearFETerm(a,itrian,iquad,[:right,:left,:left]) # We indicate here in which part of the boundary we want to define each field

This should be pretty easy to implement. @santiagobadia what do you think about this new API?

Another enhancement would be to try to merge all these terms inside the same kernel since now they lead to an enormous operation tree (you have already seen the insane stack trace...)

[santiagobadia]

I would prefer to stick to the API we proposed before, with inner and outer (or left and right, which I would not use because it is not name-revealing). Since the unknowns are defined in both sides, where is the problem using inner or outer to define the interface jumps? We should be able to evaluate both. Why is it different from standard DG?

We can probably have a chat, in order to undestand the problem that makes such approach not an option.

[fverdugo]

@oriolcg I have fixed the bug related with the computation of jumps on the FSI interface. It is available in Gridap#master, not yet released. With this, you should able to write

 function nitsche_Γ(x,y) 
   us_Γ, uf_Γ, p_Γ = x 
   vs_Γ, vf_Γ, q_Γ = y 
   uf = jump(uf_Γ)
   p = jump(p_Γ) 
   us = -jump(us_Γ) 
   vf = jump(vf_Γ) 
   q = jump(q_Γ)
   vs = -jump(vs_Γ)
   ε_uf = jump(ε(uf_Γ))
   ε_us = -jump(ε(us_Γ))
   ε_vf = jump(ε(vf_Γ))
   ε_vs = -jump(ε(vs_Γ))

 # + definition of terms here

end

You can use this for the moment until we have an alternative API.

[fverdugo]

I would prefer to stick to the API we proposed before, with inner and outer (or left and right, which I would not use because it is not name-revealing). Since the unknowns are defined in both sides, where is the problem using inner or outer to define the interface jumps? We should be able to evaluate both. Why is it different from standard DG?

The problem is that the array of cell contributions provided by nitsche_Γ(x,y) has to be complemented with other arrays that say to which cell in the underlying mesh each contribution goes to. This info is handled by the Triangulation object inside the FETerm. In particular, an SkeletonTriangulation (which is the type provided by the InterfaceTriangulation constructor) stores a pair with the ids of the cells that are on the left and right of the interface. Thus, the nitsche_Γ(x,y) has also to return something that represents a "pair", which is what happens when you implement the function based on jump, but not if you take left/right inside the function.

Implementing something like

function a(x,y)
# these would be  already in the left / rigth side
us,uf,p = x
vs,vf,q = y
# your terms here
end
t = LinearFETerm(a,itrian,iquad,[:right,:left,:left])

Is straight-forward with the current framework and the definition of a is much more compact. You only need to implement a new FETerm that internally takes left/right using the info provided in the last argument.

Implementing something like

function a(x,y)
   us_Γ, uf_Γ, p_Γ = x 
   vs_Γ, vf_Γ, q_Γ = y 
   uf = uf_Γ.left
   p = p_Γ.left 
   us = us_Γ.right 
   vf = vf_Γ.left 
   q = q_Γ.left 
   vs = vs_Γ.right
# your terms here
end
t = LinearFETerm(a,itrian,iquad)

is orders of magnitude more complicated and the definition of a is much more bloated.

That's why I was suggesting option 1. What is your concern with option 1?

[oriolcg]

@fverdugo I implemented the temporary solution you proposed with the jumps. It didn't work directly, because the operation sum (+) is not supported by the jump operator. However, the following work-around worked:

function nitsche_Γ(x,y)
  us_Γ, uf_Γ, p_Γ = x
  vs_Γ, vf_Γ, q_Γ = y
  uf = jump(uf_Γ)
  p = jump(p_Γ)
  us = -jump(us_Γ)
  vf = jump(vf_Γ)
  q = jump(q_Γ)
  vs = -jump(vs_Γ)
  εuf = 0.5 * ( jump(∇(uf_Γ)) + jump(transpose(∇(uf_Γ))) )
  εvf = 0.5 * ( jump(∇(vf_Γ)) + jump(transpose(∇(vf_Γ))) )
  εus = 0.5 * ( -jump(∇(us_Γ)) + -jump(transpose(∇(us_Γ))) )
  εvs = 0.5 * ( -jump(∇(vs_Γ)) + -jump(transpose(∇(vs_Γ))) )
  # ...
end

I also replaced the model by the widely used elastic flag after a cylinder, generated in GMSH.

I still need to add the computation of surface forces. After that, I think that the steady version of the tutorial will be ready.

[fverdugo]

@oriolcg Very nice problem setup!

@fverdugo I implemented the temporary solution you proposed with the jumps. It didn't work directly, because the operation sum (+) is not supported by the jump operator. However, the following work-around worked:

In theory, you should be able to use the symmetric gradient operator inside the jump function, i.e. jump(ε(uf_Γ)) just like the standard gradient. If not, it is a bug that we need to fix. Can you provide a small reproducer that triggers the errors you have found?

[oriolcg]

@fverdugo I added the test/fsi_failing_test.jl to reproduce the error I was observing. As it is now it's failing, but uncommenting the following lines works: https://github.com/gridap/Tutorials/blob/c0cc7f5ae2d4d3cb7f75bdc8f03573445e2c7935/test/fsi_failing_test.jl#L96-L97

[fverdugo]

@oriolcg I have fixed the bug associated with the symmetric gradient in Gridap#master. Now, you should be able to write εuf = jump(ε(uf_Γ)).

BTW, I have seen in the Manifest.toml file that you have added Gridap in development mode. This is problematic since this manifest file only works for your local machine. I suggest to track the master branch of Gridap instead:

(Tutorials) pkg> add Gridap#master

Then, we can also use the same manifest file to make sure that we are using the same versions of all dependencies.

Another remark. The master branch in the Tutorials repo has changed quite a lot since we created the fsi_tutorial branch. Could you do git pull origin master?

[oriolcg]

I just pushed an "almost ready to merge" version of the fsi_tutorial. I have added a couple of comments/questions for @fverdugo. If you have some time, please take a look at the tutorial and let me know what do you think and if there is something missing that should be added.

[fverdugo]

@oriolcg nice work!

• First question: we need a way of avoiding to use the jump function in this context. To think a nice API to achieve this. But for the moment we can use jump without detailing why we use the jump since it will change. In addition, the writing of the terms associated with the interface transmission conditions should also be simplified. See for instance, what I have achieved for embedded interfaces here https://github.com/gridap/GridapEmbedded.jl/blob/2e4b7eda0c37c14858ee12921fa3d5be594d94e6/test/GridapEmbeddedTests/BimaterialPoissonCutFEMTests.jl#L146
• Second question: I have implemented cell_measure, but not cell_size. If you want, we can have a zoom talk and comment how you can implement it (should be easy).

Some comments:

• Have you built the webpage of the tutorial locally? There are several displaying issues is the web version.

[oriolcg]

First question: we need a way of avoiding to use the jump function in this context. To think a nice API to achieve this. But for the moment we can use jump without detailing why we use the jump since it will change. In addition, the writing of the terms associated with the interface transmission conditions should also be simplified. See for instance, what I have achieved for embedded interfaces here https://github.com/gridap/GridapEmbedded.jl/blob/2e4b7eda0c37c14858ee12921fa3d5be594d94e6/test/GridapEmbeddedTests/BimaterialPoissonCutFEMTests.jl#L146

Ok. I removed the comment. The way you have for defining this term in your example looks much better than what we have in this tutorial. We cannot do the same here, right?

Have you built the webpage of the tutorial locally? There are several displaying issues is the web version.

Yes. Is there any other way of doing it?

[fverdugo]

The way you have for defining this term in your example looks much better than what we have in this tutorial. We cannot do the same here, right?

I am thinking of a way of achieving a similar API for body-fitted interfaces. It is possible, but would require some significant work in Gridap.

This is an example of the API I am thinking of. The same term as in the link above ported to the body fitted case, would be as follows (API not yet working)

function jump_u(u)
  u1,u2 = u
  u1.inwards - u2.outwards
end

function mean_q(u)
  u1,u2 = u
  κ1*α1*n_Γ*∇(u1.inwards) + κ2*α2*n_Γ*∇(u2.outwards)
end

function mean_u(u)
  u1,u2 = u
  κ2*u1.inwards + κ1*u2.outwards
end
function a_Γ(u,v)
  β*jump_u(v)*jump_u(u) - jump_u(v)*mean_q(u) - mean_q(v)*jump_u(u)
end

function l_Γ(v)
  mean_u(v)*(n_Γ*j)
end

t_Γ = AffineFETerm(a_Γ,l_Γ,trian_Γ,quad_Γ)

where in this case trian_Γ is an instance of SkeletonTriangulation (constructed via the InterfaceTriangulation constructor)

[fverdugo]

closed via PR #19