gridap / Gridap.jl

Grid-based approximation of partial differential equations in Julia
MIT License
695 stars 95 forks source link

Applying Gridap to solving problems in solid mechanics #318

Closed bhaveshshrimali closed 4 years ago

bhaveshshrimali commented 4 years ago

Hi @santiagobadia @fverdugo ,

Firstly, thank you for bringing Gridap to fruition. I am looking forward to pick it up soon. :) I am particularly targetting applications involving PDEs (and later on coupled- PDE-ODEs) with periodic boundary conditions. These occur frequently in numerical homogenization of elliptic operators.

I have a couple of questions:

I would also be happy to contribute tutorials once I make some progress on it.

Thanks again

santiagobadia commented 4 years ago

Hi @bhaveshshrimali

We have big plans for Gridap in computational solid mechanics. You have probably checked the solid mechanics tutorials. We have one example with state-depedent variables, in order to check that our machinery works for these problems. And we would like to explore a GridapFE2 module, re-implement here our octree-based additive manufacturing solvers, etc. It would be nice to have more people involved, and even better application-experts collaborators.

Gridap already supports periodic boundary conditions. See this MHD driver. It is as simple as

model = CartesianDiscreteModel(domain,partition; isperiodic=(false,false,true))

where (false,false,true) means periodic in Z direction and not in X and Y, as you can guess.

We already have the concept of (multi-dof) linear constraints, you can create a FESpaceWithLinearConstraints providing these constraints in a particular way. The Gridap assembler can deal with them. We plan to use this soon for h-adaptivity in finite elements. But I think that with the previous method you don't need to go that deep.

Thanks for the interest.

fverdugo commented 4 years ago

@bhaveshshrimali

Here you have a minimal working example with periodic BCs (Poisson Eq.)

using Gridap

u(x) = sin(2*pi*(x[1]-0.25))*x[2]
f(x) = -Δ(u)(x)

domain = (0,1,0,2)
cells = (10,20)
model = CartesianDiscreteModel(domain,cells;isperiodic=(true,false))

labels = get_face_labeling(model)
add_tag_from_tags!(labels,"dirichlet",[1,2,3,4,5,6])

order = 2
V = TestFESpace(
 reffe=:Lagrangian, order=order, valuetype=Float64,
 conformity=:H1, model=model, dirichlet_tags="dirichlet")

U = TrialFESpace(V,u)

trian = Triangulation(model)
quad = CellQuadrature(trian,2*order)

a(u,v) = ∇(v)⋅∇(u)
l(v) = v*f
t_Ω = AffineFETerm(a,l,trian,quad)

op = AffineFEOperator(U,V,t_Ω)
uh = solve(op)

writevtk(trian,"pcbs",order=order,cellfields=["uh"=>uh])

pbcs

Hope it helps.

bhaveshshrimali commented 4 years ago

Thanks very much to both of you.

You have probably checked the solid mechanics tutorials. We have one example with state-depedent variables, in order to check that our machinery works for these problems.

Thanks @santiagobadia . I started with the Hyperelasticity example and the [damage example] (https://gridap.github.io/Tutorials/stable/pages/t010_isotropic_damage/) is my next on my list.

Here you have a minimal working example with periodic BCs (Poisson Eq.)

Thanks @fverdugo This is perfect

Appreciate both your responses. Will keep up with the latest developments on this. Thanks again!

avigliotti commented 4 years ago

if you guys have interest in solid mechanics you might also want to give a look at this

https://github.com/avigliotti/AD4SM.jl

it uses automatic differentiation to state and solve solid mechanics problems, solution is found as the minima of free energy, with residual forces being the gradient and tangent stiffness matrix the hessian of structure's strain energy, respectively.

p.s. didn't mean to steal focus from Gridap.jl, just share some work

santiagobadia commented 4 years ago

Thanks @avigliotti for sharing this info.

We have been able to compute jacobians out of residuals and also jacobians and residuals out of energy functionals using ForwardDiff dual-number bassed AD engine. But for much simpler models than the ones you have considered.

Hopefully, we will create tutorials with these capabilities soon.

avigliotti commented 4 years ago

The implementation of automatic differentiation in AD4SM is also based on ForwardDiff, which is a terrific and elegant package, the difference is that in AD4SM the Hessian is not evaluated recursively as the gradient of the gradient's components, but it is treated explicitly to take advantage of Hessian's symmetry, this way we only compute the upper half of it

santiagobadia commented 4 years ago

Thanks for all your comments.

I will close the issue, but feel free to open it again or create another one if needed.