-
Notifications
You must be signed in to change notification settings - Fork 37
/
search_index.js
3 lines (3 loc) · 197 KB
/
search_index.js
1
2
3
var documenterSearchIndex = {"docs":
[{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"EditURL = \"https://github.com/gridap/Tutorials/blob/master/src/validation.jl\"","category":"page"},{"location":"pages/t002_validation/#Tutorial-2:-Code-validation-1","page":"2 Code validation","title":"Tutorial 2: Code validation","text":"","category":"section"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"(Image: ) (Image: )","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"In this tutorial, we will learn","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"How to implement the method of manufactured solutions\nHow to perform a convergence test\nHow to define the discretization error\nHow to integrate error norms\nHow to generate Cartesian meshes in arbitrary dimensions","category":"page"},{"location":"pages/t002_validation/#Problem-statement-1","page":"2 Code validation","title":"Problem statement","text":"","category":"section"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"In this tutorial, we show how to validate a code using the well known method of manufactured solutions. For the sake of simplicity, we consider the Poisson equation in the unit square Omegadoteq (01)^2 as a model problem,","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"leftlbrace\nbeginaligned\n-Delta u = f textin Omega\nu = g texton partialOmega\nendaligned\nright","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"We are going to consider two different manufactured solutions. On the one hand, we consider function u(x)=x_1+x_2, which can be exactly represented by the FE interpolation that we construct below. Thus, one expects that the obtained approximation error is near the machine precision. We are going to check that this is true in the code. On the other hand, we consider a function that cannot be captured exactly by the interpolation, namely u(x)=x_2 sin(2 pi x_1). Here, our goal is to confirm that the convergence order of the discretization error is the optimal one.","category":"page"},{"location":"pages/t002_validation/#Manufactured-solution-1","page":"2 Code validation","title":"Manufactured solution","text":"","category":"section"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"We start by defining the manufactured solution u(x) = x_1+x_2 and the source term f associated with it, namely fdoteq-Delta(x_1+x_2)=0.","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"using Gridap\n\nu(x) = x[1] + x[2]\nf(x) = 0","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"We also need to define the gradient of u since we will compute the H^1 error norm later. In that case, the gradient is simply defined as","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"∇u(x) = VectorValue(1,1)","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"Note that we have used the constructor VectorValue to build the vector that represents the gradient. However, we still need a final trick. We need to tell the Gridap library that the gradient of the function u is available in the function ∇u (at this moment u and ∇u are two standard Julia functions without any connection between them). This is done by adding an extra method to the function gradient (aka ∇) defined in Gridap:","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"import Gridap: ∇\n∇(::typeof(u)) = ∇u","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"Now, it is possible to recover function ∇u from function u as ∇(u). You can check that the following expression evaluates to true.","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"∇(u) === ∇u","category":"page"},{"location":"pages/t002_validation/#Cartesian-mesh-generation-1","page":"2 Code validation","title":"Cartesian mesh generation","text":"","category":"section"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"In order to discretize the geometry of the unit square, we use the Cartesian mesh generator available in Gridap:","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"domain = (0,1,0,1)\npartition = (4,4)\nmodel = CartesianDiscreteModel(domain,partition)","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"The type CartesianDiscreteModel is a concrete type that inherits from DiscreteModel, which is specifically designed for building Cartesian meshes. The CartesianDiscreteModel constructor takes a tuple containing limits of the box we want to discretize plus a tuple with the number of cells to be generated in each direction (here 4 by 4 cells). Note that the CaresianDiscreteModel is implemented for arbitrary dimensions. For instance, the following lines build a CartesianDiscreteModel for the unit cube (01)^3 with 4 cells per direction","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"domain3d = (0,1,0,1,0,1)\npartition3d = (4,4,4)\nmodel3d = CartesianDiscreteModel(domain3d,partition3d)","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"You could also generate a mesh for the unit tesseract (01)^4 (i.e., the unit cube in 4D). Look how the 2D and 3D models are built and just follow the sequence.","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"Let us return to the 2D CartesianDiscreteModel that we have already constructed. You can inspect it by writing it into vtk format. Note that you can also print a 3D model, but not a 4D one. In the future, it would be cool to generate a movie from a 4D model, but this functionality is not yet implemented.","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"writevtk(model,\"model\")","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"If you open the generated files, you will see that the boundary vertices and facets are identified with the name \"boundary\". This is just what we need to impose the Dirichlet boundary conditions in this example.","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"These are the vertices in the model","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"(Image: )","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"and these the facets","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"(Image: )","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"As you can see in the pictures, the objects on the boundary are correctly tagged with the name \"boundary\".","category":"page"},{"location":"pages/t002_validation/#FE-approximation-1","page":"2 Code validation","title":"FE approximation","text":"","category":"section"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"We compute a FE approximation of the Poisson problem above by following the steps detailed in the previous tutorial:","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"order = 1\nV0 = TestFESpace(\n reffe=:Lagrangian, order=order, valuetype=Float64,\n conformity=:H1, model=model, dirichlet_tags=\"boundary\")\n\nU = TrialFESpace(V0,u)\n\ntrian = Triangulation(model)\ndegree = 2\nquad = CellQuadrature(trian,degree)\n\na(u,v) = ∇(v)⊙∇(u)\nb(v) = v*f\n\nt_Ω = AffineFETerm(a,b,trian,quad)\nop = AffineFEOperator(U,V0,t_Ω)\n\nuh = solve(op)","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"Note that we are imposing Dirichlet boundary conditions on the objects tagged as \"boundary\" and that we are using the manufactured solution u to construct the trial FE space. Not also that we are not explicitly constructing an Assembler object nor a FESolver. We are relying on default values.","category":"page"},{"location":"pages/t002_validation/#Measuring-the-discretization-error-1","page":"2 Code validation","title":"Measuring the discretization error","text":"","category":"section"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"Our goal is to check that the discratization error associated with the computed approximation uh is close to machine precision. To this end, the first step is to compute the discretization error, which is done as you would expect:","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"e = u - uh","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"Once the error is defined, you can, e.g., visualize it.","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"writevtk(trian,\"error\",cellfields=[\"e\" => e])","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"This generates a file called error.vtu. Open it with Paraview to check that the error is of the order of the machine precision.","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"(Image: )","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"A more rigorous way of quantifying the error is to measure it with a norm. Here, we use the L^2 and H^1 norms, which are defined as","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":" w _L^2^2 doteq int_Omega w^2 textdOmega quad\n w _H^1^2 doteq int_Omega w^2 + nabla w cdot nabla w textdOmega","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"In order to compute these norms, we are going to use the integrate function. To this end, we need to define the integrands that we want to integrate, namely","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"l2(w) = w*w\nh1(w) = a(w,w) + l2(w)","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"Note that in order to define the integrand of the H^1 norm, we have reused function a, previously used to define the bilinear form of the problem. Once we have defined the integrands, we are ready to compute the integrals. For the L^2 norm","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"el2 = sqrt(sum( integrate(l2(e),trian,quad) ))","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"and for the H^1 norm","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"The integrate function works as follows. In the first argument, we pass the integrand. In the second and third arguments, we pass a Triangulation object and aCellQuadrature that represent the data needed in order to perform the integrals numerically. The integrate function returns an object containing the contribution to the integrated value of each cell in the given Triangulation. To end up with the desired error norms, one has to sum these contributions and take the square root. You can check that the computed error norms are close to machine precision (as one would expect).","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"tol = 1.e-10\n@assert el2 < tol\n@assert eh1 < tol","category":"page"},{"location":"pages/t002_validation/#Convergence-test-1","page":"2 Code validation","title":"Convergence test","text":"","category":"section"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"We end up this tutorial by performing a convergence test, where we are going to use all the new concepts we have learned. We will consider a manufactured solution that does not belong to the FE interpolation space. In this test, we expect to see the optimal convergence order of the FE discretization.","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"Here, we define the manufactured functions","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"const k = 2*pi\nu(x) = sin(k*x[1]) * x[2]\n∇u(x) = VectorValue(k*cos(k*x[1])*x[2], sin(k*x[1]))\nf(x) = (k^2)*sin(k*x[1])*x[2]","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"Since we have redefined the valiables u, ∇u, and f, we need to execute these lines again","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"∇(::typeof(u)) = ∇u\nb(v) = v*f","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"In order to perform the convergence test, we write in a function all the code needed to perform a single computation and measure its error. The input of this function is the number of cells in each direction and the interpolation order. The output is the computed L^2 and H^1 error norms.","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"function run(n,order)\n\n domain = (0,1,0,1)\n partition = (n,n)\n model = CartesianDiscreteModel(domain,partition)\n\n V0 = TestFESpace(\n reffe=:Lagrangian, order=order, valuetype=Float64,\n conformity=:H1, model=model, dirichlet_tags=\"boundary\")\n\n U = TrialFESpace(V0,u)\n\n trian = Triangulation(model)\n degree=2*order\n quad = CellQuadrature(trian,degree)\n\n t_Ω = AffineFETerm(a,b,trian,quad)\n op = AffineFEOperator(U,V0,t_Ω)\n\n uh = solve(op)\n\n e = u - uh\n\n el2 = sqrt(sum( integrate(l2(e),trian,quad) ))\n eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))\n\n (el2, eh1)\n\nend","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"The following function does the convergence test. It takes a vector of integers (representing the number of cells per direction in each computation) plus the interpolation order. It returns the L^2 and H^1 error norms for each computation as well as the corresponding cell size.","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"function conv_test(ns,order)\n\n el2s = Float64[]\n eh1s = Float64[]\n hs = Float64[]\n\n for n in ns\n\n el2, eh1 = run(n,order)\n h = 1.0/n\n\n push!(el2s,el2)\n push!(eh1s,eh1)\n push!(hs,h)\n\n end\n\n (el2s, eh1s, hs)\n\nend","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"We are ready to perform the test! We consider several mesh sizes and interpolation order equal to 2.","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"el2s, eh1s, hs = conv_test([8,16,32,64,128],2);\nnothing #hide","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"With the generated data, we do the classical convergence plot.","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"using Plots\n\nplot(hs,[el2s eh1s],\n xaxis=:log, yaxis=:log,\n label=[\"L2\" \"H1\"],\n shape=:auto,\n xlabel=\"h\",ylabel=\"error norm\")","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"If you run the code in a notebook, you will see a figure like this one: (Image: )","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"The generated curves make sense. It is observed that the convergence of the H^1 error is slower that L^2 one. However, in order to be more conclusive, we need to compute the slope of these lines. It can be done with this little function that internally uses a linear regression.","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"function slope(hs,errors)\n x = log10.(hs)\n y = log10.(errors)\n linreg = hcat(fill!(similar(x), 1), x) \\ y\n linreg[2]\nend","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"The slope for the L^2 error norm is computed as","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"slope(hs,el2s)","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"and for the H^1 error norm","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"slope(hs,eh1s)","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"If your run these lines in a notebook, you will see that the slopes for the L^2 and H^1 error norms are circa 3 and 2 respectively (as one expects for interpolation order 2)","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"Congrats, another tutorial done!","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"","category":"page"},{"location":"pages/t002_validation/#","page":"2 Code validation","title":"2 Code validation","text":"This page was generated using Literate.jl.","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"EditURL = \"https://github.com/gridap/Tutorials/blob/master/src/hyperelasticity.jl\"","category":"page"},{"location":"pages/t005_hyperelasticity/#Tutorial-5:-Hyper-elasticity-1","page":"5 Hyper-elasticity","title":"Tutorial 5: Hyper-elasticity","text":"","category":"section"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"(Image: ) (Image: )","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"note: Note\nThis tutorial is under construction, but the code below is already functional.","category":"page"},{"location":"pages/t005_hyperelasticity/#Problem-statement-1","page":"5 Hyper-elasticity","title":"Problem statement","text":"","category":"section"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"using Gridap\nusing LinearAlgebra: inv, det\nusing LineSearches: BackTracking","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Material parameters","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"const λ = 100.0\nconst μ = 1.0","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Deformation Gradient","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"F(∇u) = one(∇u) + ∇u'\n\nJ(F) = sqrt(det(C(F)))\n\n#Green strain\n\n#E(F) = 0.5*( F'*F - one(F) )\n\n@law dE(∇du,∇u) = 0.5*( ∇du⋅F(∇u) + (∇du⋅F(∇u))' )","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Right Cauchy-green deformation tensor","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"C(F) = (F')⋅F","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Constitutive law (Neo hookean)","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"@law function S(∇u)\n Cinv = inv(C(F(∇u)))\n μ*(one(∇u)-Cinv) + λ*log(J(F(∇u)))*Cinv\nend\n\n@law function dS(∇du,∇u)\n Cinv = inv(C(F(∇u)))\n _dE = dE(∇du,∇u)\n\tλ*(Cinv⊙_dE)*Cinv + 2*(μ-λ*log(J(F(∇u))))*Cinv⋅_dE⋅(Cinv')\nend","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Cauchy stress tensor","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"@law σ(∇u) = (1.0/J(F(∇u)))*F(∇u)⋅S(∇u)⋅(F(∇u))'","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Weak form","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"res(u,v) = dE(∇(v),∇(u)) ⊙ S(∇(u))\n\njac_mat(u,du,v) = dE(∇(v),∇(u)) ⊙ dS(∇(du),∇(u))\n\njac_geo(u,du,v) = ∇(v) ⊙ ( S(∇(u))⋅∇(du) )\n\njac(u,du,v) = jac_mat(u,v,du) + jac_geo(u,v,du)","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Model","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"domain = (0,1,0,1)\npartition = (20,20)\nmodel = CartesianDiscreteModel(domain,partition)","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Define new boundaries","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"labels = get_face_labeling(model)\nadd_tag_from_tags!(labels,\"diri_0\",[1,3,7])\nadd_tag_from_tags!(labels,\"diri_1\",[2,4,8])","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Construct the FEspace","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"V = TestFESpace(\n model=model,valuetype=VectorValue{2,Float64},\n reffe=:Lagrangian,conformity=:H1,order=1,\n dirichlet_tags = [\"diri_0\", \"diri_1\"])","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Setup integration","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"trian = Triangulation(model)\ndegree = 2\nquad = CellQuadrature(trian,degree)","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Setup weak form terms","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"t_Ω = FETerm(res,jac,trian,quad)","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Setup non-linear solver","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"nls = NLSolver(\n show_trace=true,\n method=:newton,\n linesearch=BackTracking())\n\nsolver = FESolver(nls)\n\nfunction run(x0,disp_x,step,nsteps)\n\n g0 = VectorValue(0.0,0.0)\n g1 = VectorValue(disp_x,0.0)\n U = TrialFESpace(V,[g0,g1])\n\n #FE problem\n op = FEOperator(U,V,t_Ω)\n\n println(\"\\n+++ Solving for disp_x $disp_x in step $step of $nsteps +++\\n\")\n\n uh = FEFunction(U,x0)\n\n uh, = solve!(uh,solver,op)\n\n writevtk(trian,\"results_$(lpad(step,3,'0'))\",cellfields=[\"uh\"=>uh,\"sigma\"=>σ(∇(uh))])\n\n return get_free_values(uh)\n\nend\n\nfunction runs()\n\n disp_max = 0.75\n disp_inc = 0.02\n nsteps = ceil(Int,abs(disp_max)/disp_inc)\n\n x0 = zeros(Float64,num_free_dofs(V))\n\n for step in 1:nsteps\n disp_x = step * disp_max / nsteps\n x0 = run(x0,disp_x,step,nsteps)\n end\n\nend\n\n#Do the work!\nruns()","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Picture of the last load step (Image: )","category":"page"},{"location":"pages/t005_hyperelasticity/#Extension-to-3D-1","page":"5 Hyper-elasticity","title":"Extension to 3D","text":"","category":"section"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"Extending this tutorial to the 3D case is straightforward. It is leaved as an exercise.","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"(Image: )","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"","category":"page"},{"location":"pages/t005_hyperelasticity/#","page":"5 Hyper-elasticity","title":"5 Hyper-elasticity","text":"This page was generated using Literate.jl.","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"EditURL = \"https://github.com/gridap/Tutorials/blob/master/src/darcy.jl\"","category":"page"},{"location":"pages/t007_darcy/#Tutorial-7:-Darcy-equation-(with-RT)-1","page":"7 Darcy equation (with RT)","title":"Tutorial 7: Darcy equation (with RT)","text":"","category":"section"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"(Image: ) (Image: )","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"In this tutorial, we will learn","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"How to implement multi-field PDEs\nHow to build div-conforming FE spaces\nHow to impose boundary conditions in multi-field problems","category":"page"},{"location":"pages/t007_darcy/#Problem-statement-1","page":"7 Darcy equation (with RT)","title":"Problem statement","text":"","category":"section"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"In this tutorial, we show how to solve a multi-field PDE in Gridap. As a model problem, we consider the Darcy equations with Dirichlet and Neumann boundary conditions. The PDE we want to solve is: find the flux vector u, and the pressure p such that","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":" leftlbrace\n beginaligned\n kappa^-1 u + nabla p = 0 textin Omega\n nabla cdot u = f textin Omega\n u cdot n = g texton Gamma_rm D\n p = h texton Gamma_rm N\n endaligned\n right","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"being n the outwards unit normal vector to the boundary partialOmega. In this particular tutorial, we consider the unit square Omega doteq (01)^2 as the computational domain, the Neumann boundary Gamma_rm N is the right and left sides of Omega, and Gamma_rm D is the bottom and top sides of Omega. We consider f = g doteq 0 and h(x) doteq x_1, i.e., h equal to 0 on the left side and 1 on the right side. The inverse of the permeability tensor, namely kappa^-1(x), is chosen equal to","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"beginpmatrix\n 100 90 \n 90 100\nendpmatrix\ntext for x in 0406^2 text and \nbeginpmatrix\n 1 0 \n 0 1\nendpmatrix\n text\totherwise","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"In order to state this problem in weak form, we introduce the following Sobolev spaces. H(mathrmdivOmega) is the space of vector fields in Omega, whose components and divergence are in L^2(Omega). On the other hand, H_g(mathrmdivOmega) and H_0(mathrmdivOmega) are the subspaces of functions in H(mathrmdivOmega) such that their normal traces are equal to g and 0 respectively almost everywhere in Gamma_rm D. With these notations, the weak form reads: find (up)in H_g(mathrmdivOmega)times L^2(Omega) such that a((uq)(vq)) = b(vq) for all (vq)in H_0(mathrmdivOmega)times L^2(Omega), where","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"beginaligned\na((uq)(vq)) doteq int_Omega v cdot left(kappa^-1 uright) rm dOmega - int_Omega (nabla cdot v) p rm dOmega + int_Omega q (nabla cdot u) rm dOmega\nb(vq) doteq int_Omega q f rm dOmega - int_Gamma_rm N (vcdot n) h rm dGamma\nendaligned","category":"page"},{"location":"pages/t007_darcy/#Numerical-scheme-1","page":"7 Darcy equation (with RT)","title":"Numerical scheme","text":"","category":"section"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"In this tutorial, we use the Raviart-Thomas (RT) space for the flux approximation [1]. On a reference square with sides aligned with the Cartesian axes, the \\ac{rt} space of order k is represented as Q_(k+1k) times Q_(kk+1), being the polynomial space defined as follows. The component w_alpha of a vector field w in Q_(k+1k) times Q_(kk+1) is obtained as the tensor product of univariate polynomials of order k+1 in direction alpha times univariate polynomials of order k on the other directions. That is, nablacdot w in Q_k, where Q_k is the multivariate polynomial space of degree at most k in each of the spatial coordinates. Note that the definition of the \\ac{rt} space also applies to arbitrary dimensions. The global FE space for the flux V is obtained by mapping the cell-wise RT space into the physical space using the Piola transformation and enforcing continuity of normal traces across cells (see [1] for specific details).","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"We consider the subspace V_0 of functions in V with zero normal trace on Gamma_rm D, and the subspace V_g of functions in V with normal trace equal to the projection of g onto the space of traces of V on Gamma_rm D. With regard to the pressure, we consider the discontinuous space of cell-wise polynomials in Q_k.","category":"page"},{"location":"pages/t007_darcy/#Discrete-model-1","page":"7 Darcy equation (with RT)","title":"Discrete model","text":"","category":"section"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"We start the driver loading the Gridap package and constructing the geometrical model. We generate a 100times100 structured mesh for the domain (01)^2.","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"using Gridap\ndomain = (0,1,0,1)\npartition = (100,100)\nmodel = CartesianDiscreteModel(domain,partition)","category":"page"},{"location":"pages/t007_darcy/#Multi-field-FE-spaces-1","page":"7 Darcy equation (with RT)","title":"Multi-field FE spaces","text":"","category":"section"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"Next, we build the FE spaces. We consider the first order RT space for the flux and the discontinuous pressure space as described above. This mixed FE pair satisfies the inf-sup condition and, thus, it is stable.","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"order = 1\n\nV = TestFESpace(\n reffe=:RaviartThomas, order=order, valuetype=VectorValue{2,Float64},\n conformity=:HDiv, model=model, dirichlet_tags=[5,6])\n\nQ = TestFESpace(\n reffe=:QLagrangian, order=order, valuetype=Float64,\n conformity=:L2, model=model)","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"Note that the Dirichlet boundary for the flux are the bottom and top sides of the squared domain (identified with the boundary tags 5, and 6 respectively), whereas no Dirichlet data can be imposed on the pressure space. We select conformity=:HDiv for the flux (i.e., shape functions with H^1(mathrmdivOmega) regularity) and conformity=:L2 for the pressure (i.e. discontinuous shape functions).","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"From these objects, we construct the trial spaces. Note that we impose homogeneous boundary conditions for the flux.","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"uD = VectorValue(0.0,0.0)\nU = TrialFESpace(V,uD)\nP = TrialFESpace(Q)","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"When the singe-field spaces have been designed, the multi-field test and trial spaces are expressed as arrays of single-field ones in a natural way.","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"Y = MultiFieldFESpace([V, Q])\nX = MultiFieldFESpace([U, P])","category":"page"},{"location":"pages/t007_darcy/#Numerical-integration-1","page":"7 Darcy equation (with RT)","title":"Numerical integration","text":"","category":"section"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"In this example we need to integrate in the interior of Omega and on the Neumann boundary Gamma_rm N. For the volume integrals, we extract the triangulation from the geometrical model and define the corresponding cell-wise quadrature of degree of exactness at least 2 as follows.","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"trian = Triangulation(model)\ndegree = 2\nquad = CellQuadrature(trian,degree)","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"In order to integrate the Neumann boundary condition, we only need to build an integration mesh for the right side of the domain (which is the only part of Gamma_rm N, where the Neumann function h is different from zero). Within the model, the right side of Omega is identified with the boundary tag 8. Using this identifier, we extract the corresponding surface triangulation and create a quadrature with the desired degree of exactness.","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"neumanntags = [8,]\nbtrian = BoundaryTriangulation(model,neumanntags)\nbquad = CellQuadrature(btrian,degree)","category":"page"},{"location":"pages/t007_darcy/#Weak-form-1","page":"7 Darcy equation (with RT)","title":"Weak form","text":"","category":"section"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"We start by defining the permeability tensors inverses commented above using the @law macro.","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"const kinv1 = TensorValue(1.0,0.0,0.0,1.0)\nconst kinv2 = TensorValue(100.0,90.0,90.0,100.0)\n@law function σ(x,u)\n if ((abs(x[1]-0.5) <= 0.1) && (abs(x[2]-0.5) <= 0.1))\n return kinv2⋅u\n else\n return kinv1⋅u\n end\nend","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"With this definition, we can express the integrand of the bilinear form as follows.","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"px = get_physical_coordinate(trian)\n\nfunction a(x,y)\n v, q = y\n u, p = x\n v⋅σ(px,u) - (∇⋅v)*p + q*(∇⋅u)\nend","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"The arguments x and y of previous function represent a trial and a test function in the multi-field test and trial spaces X and Y respectively. In the first lines in the function definition, we unpack the single-field test and trial functions from the multi-field ones. E.g., v represents a test function for the flux and q for the pressure. These quantities can also be written as y[1] and y[2] respectively. From the single-field functions, we write the different terms of the bilinear form as we have done in previous tutorials.","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"In a similar way, we can define the forcing term related to the Neumann boundary condition.","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"nb = get_normal_vector(btrian)\nh = -1.0\nfunction b_ΓN(y)\n v, q = y\n (v⋅nb)*h\nend","category":"page"},{"location":"pages/t007_darcy/#Multi-field-FE-problem-1","page":"7 Darcy equation (with RT)","title":"Multi-field FE problem","text":"","category":"section"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"Finally, we can assemble the FE problem and solve it. Note that we build the AffineFEOperator object using the multi-field trial and test spaces Y and X.","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"t_Ω = LinearFETerm(a,trian,quad)\nt_ΓN = FESource(b_ΓN,btrian,bquad)\nop = AffineFEOperator(X,Y,t_Ω,t_ΓN)\nxh = solve(op)\nuh, ph = xh","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"Since this is a multi-field example, the solve function returns a multi-field solution xh, which can be unpacked in order to finally recover each field of the problem. The resulting single-field objects can be visualized as in previous tutorials (see next figure).","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"writevtk(trian,\"darcyresults\",cellfields=[\"uh\"=>uh,\"ph\"=>ph])","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"(Image: )","category":"page"},{"location":"pages/t007_darcy/#References-1","page":"7 Darcy equation (with RT)","title":"References","text":"","category":"section"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"[1] F. Brezzi and M. Fortin. Mixed and hybrid finite element methods. Springer-Verlag, 1991.","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"","category":"page"},{"location":"pages/t007_darcy/#","page":"7 Darcy equation (with RT)","title":"7 Darcy equation (with RT)","text":"This page was generated using Literate.jl.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"EditURL = \"https://github.com/gridap/Tutorials/blob/master/src/dg_discretization.jl\"","category":"page"},{"location":"pages/t006_dg_discretization/#Tutorial-6:-Poisson-equation-(with-DG)-1","page":"6 Poisson equation (with DG)","title":"Tutorial 6: Poisson equation (with DG)","text":"","category":"section"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"(Image: ) (Image: )","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"In this tutorial, we will learn","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"How to solve a simple PDE with a DG method\nHow to compute jumps and averages of quantities on the mesh skeleton\nHow to implement the method of manufactured solutions\nHow to integrate error norms\nHow to generate Cartesian meshes in arbitrary dimensions","category":"page"},{"location":"pages/t006_dg_discretization/#Problem-statement-1","page":"6 Poisson equation (with DG)","title":"Problem statement","text":"","category":"section"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"The goal of this tutorial is to solve a PDE using a DG formulation. For simplicity, we take the Poisson equation on the unit cube Omega doteq (01)^3 as the model problem, namely","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"leftlbrace\nbeginaligned\n-Delta u = f textin Omega\nu = g texton partialOmega\nendaligned\nright","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"where f is the source term and g is the prescribed Dirichlet boundary function. In this tutorial, we follow the method of manufactured solutions since we want to illustrate how to compute discretization errors. We take u(x) = 3 x_1 + x_2 + 2 x_3 as the exact solution of the problem, for which f=0 and g(x) = u(x). The selected manufactured solution u is a first order multi-variate polynomial, which can be represented exactly by the FE interpolation that we are going to define below. In this scenario, the discretization error has to be close to the machine precision. We will use this result to validate the proposed implementation.","category":"page"},{"location":"pages/t006_dg_discretization/#Numerical-Scheme-1","page":"6 Poisson equation (with DG)","title":"Numerical Scheme","text":"","category":"section"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"We consider a DG formulation to approximate the problem. In particular, we consider the symmetric interior penalty method (see, e.g. [1], for specific details). For this formulation, the approximation space is made of discontinuous piece-wise polynomials, namely","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"V doteq vin L^2(Omega) v_Tin Q_p(T) text for all TinmathcalT ","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"where mathcalT is the set of all cells T of the FE mesh, and Q_p(T) is a polynomial space of degree p defined on a generic cell T. For simplicity, we consider Cartesian meshes in this tutorial. In this case, the space Q_p(T) is made of multi-variate polynomials up to degree p in each spatial coordinate.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"In order to write the weak form of the problem, we need to introduce some notation. The sets of interior and boundary facets associated with the FE mesh mathcalT are denoted here as mathcalF_Gamma and mathcalF_partialOmega respectively. In addition, for a given function vin V restricted to the interior facets mathcalF_Gamma, we introduce the well known jump and mean value operators,","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"beginaligned\nlbracklbrack v n rbrackrbrack doteq v^+ n^+ + v^- n^-\n nabla v doteq dfrac nabla v^+ + nabla v^-2\nendaligned","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"with v^+, and v^- being the restrictions of vin V to the cells T^+, T^- that share a generic interior facet in mathcalF_Gamma, and n^+, and n^- are the facet outward unit normals from either the perspective of T^+ and T^- respectively.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"With this notation, the weak form associated with the interior penalty formulation of our problem reads: find uin V such that a(uv) = b(v) for all vin V. The bilinear and linear forms a(cdotcdot) and b(cdot) have contributions associated with the bulk of Omega, the boundary facets mathcalF_partialOmega, and the interior facets mathcalF_Gamma, namely","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"beginaligned\na(uv) = a_Omega(uv) + a_partialOmega(uv) + a_Gamma(uv)\nb(v) = b_Omega(v) + b_partialOmega(v)\nendaligned","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"These contributions are defined as","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"beginaligned\na_Omega(uv) doteq sum_TinmathcalT int_T nabla v cdot nabla u rm dT\nb_Omega(v) doteq int_Omega v f rm dOmega\nendaligned","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"for the volume,","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"beginaligned\na_partialOmega(uv) doteq sum_FinmathcalF_partialOmega dfracgammaF int_F v u rm dF - sum_FinmathcalF_partialOmega int_F v (nabla u cdot n) rm dF - sum_FinmathcalF_partialOmega int_F (nabla v cdot n) u rm dF \nb_partialOmega doteq sum_FinmathcalF_partialOmega dfracgammaF int_F v g rm dF - sum_FinmathcalF_partialOmega int_F (nabla v cdot n) g rm dF\nendaligned","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"for the boundary facets and,","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"beginaligned\na_Gamma(uv) doteq sum_FinmathcalF_Gamma dfracgammaF int_F lbracklbrack v n rbrackrbrackcdot lbracklbrack u n rbrackrbrack rm dF - sum_FinmathcalF_Gamma int_F lbracklbrack v n rbrackrbrackcdot nabla u rm dF - sum_FinmathcalF_Gamma int_F nabla v cdot lbracklbrack u n rbrackrbrack rm dF\nendaligned","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"for the interior facets. In previous expressions, F denotes the diameter of the face F (in our Cartesian grid, this is equivalent to the characteristic mesh size h), and gamma is a stabilization parameter that should be chosen large enough such that the bilinear form a(cdotcdot) is stable and continuous. Here, we take gamma = p (p+1) as done in the numerical experiments in reference [2].","category":"page"},{"location":"pages/t006_dg_discretization/#Manufactured-solution-1","page":"6 Poisson equation (with DG)","title":"Manufactured solution","text":"","category":"section"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"We start by loading the Gridap library and defining the manufactured solution u and the associated source term f and Dirichlet function g.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"using Gridap\nu(x) = 3*x[1] + x[2] + 2*x[3]\nf(x) = 0\ng(x) = u(x)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"We also need to define the gradient of u since we will compute the H^1 error norm later. In that case, the gradient is simply defined as","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"∇u(x) = VectorValue(3,1,2)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"In addition, we need to tell the Gridap library that the gradient of the function u is available in the function ∇u (at this moment u and ∇u are two standard Julia functions without any connection between them). This is done by adding an extra method to the function gradient (aka ∇) defined in Gridap:","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"import Gridap: ∇\n∇(::typeof(u)) = ∇u","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"Now, it is possible to recover function ∇u from function u as ∇(u). You can check that the following expression evaluates to true.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"∇(u) === ∇u","category":"page"},{"location":"pages/t006_dg_discretization/#Cartesian-mesh-generation-1","page":"6 Poisson equation (with DG)","title":"Cartesian mesh generation","text":"","category":"section"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"In order to discretize the geometry of the unit cube, we use the Cartesian mesh generator available in Gridap.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"L = 1.0\ndomain = (0.0, L, 0.0, L, 0.0, L)\nn = 4\npartition = (n,n,n)\nmodel = CartesianDiscreteModel(domain,partition)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"The type CartesianDiscreteModel is a concrete type that inherits from DiscreteModel, which is specifically designed for building Cartesian meshes. The CartesianDiscreteModel constructor takes a tuple containing limits of the box we want to discretize plus a tuple with the number of cells to be generated in each direction (here 4times4times4 cells). You can write the model in vtk format to visualize it (see next figure).","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":" writevtk(model,\"model\")","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"(Image: )","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"Note that the CaresianDiscreteModel is implemented for arbitrary dimensions. For instance, the following lines build a CartesianDiscreteModel for the unit square (01)^2 with 4 cells per direction","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"domain2D = (0.0, L, 0.0, L)\npartition2D = (n,n)\nmodel2D = CartesianDiscreteModel(domain2D,partition2D)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"You could also generate a mesh for the unit tesseract (01)^4 (i.e., the unit cube in 4D). Look how the 2D and 3D models are built and just follow the sequence.","category":"page"},{"location":"pages/t006_dg_discretization/#FE-spaces-1","page":"6 Poisson equation (with DG)","title":"FE spaces","text":"","category":"section"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"On top of the discrete model, we create the discontinuous space V as follows","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"order = 3\nV = TestFESpace(\n reffe=:Lagrangian, valuetype=Float64, order=order,\n conformity=:L2, model=model)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"We have select a Lagrangian, scalar-valued interpolation of order 3 within the cells of the discrete model. Since the cells are hexahedra, the resulting Lagrangian shape functions are tri-cubic polynomials. In contrast to previous tutorials, where we have constructed H^1-conforming (i.e., continuous) FE spaces, here we construct a L^2-conforming (i.e., discontinuous) FE space. That is, we do not impose any type of continuity of the shape function on the cell boundaries, which leads to the discontinuous FE space V of the DG formulation. Note also that we do not pass any information about the Dirichlet boundary to the TestFESpace constructor since the Dirichlet boundary conditions are not imposed strongly in this example.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"From the V object we have constructed in previous code snippet, we build the trial FE space as usual.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"U = TrialFESpace(V)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"Note that we do not pass any Dirichlet function to the TrialFESpace constructor since we do not impose Dirichlet boundary conditions strongly here.","category":"page"},{"location":"pages/t006_dg_discretization/#Numerical-integration-1","page":"6 Poisson equation (with DG)","title":"Numerical integration","text":"","category":"section"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"Once the FE spaces are ready, the next step is to set up the numerical integration. In this example, we need to integrate in three different domains: the volume covered by the cells mathcalT (i.e., the computational domain Omega), the surface covered by the boundary facets mathcalF_partialOmega (i.e., the boundary partialOmega), and the surface covered by the interior facets mathcalF_Gamma (i.e. the so-called mesh skeleton). In order to integrate in Omega and on its boundary partialOmega, we use Triangulation and BoundaryTriangulation objects as already discussed in previous tutorials.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"trian = Triangulation(model)\nbtrian = BoundaryTriangulation(model)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"Here, we do not pass any boundary identifier to the BoundaryTriangulation constructor. In this case, an integration mesh for the entire boundary partialOmega is constructed by default (which is just what we need in this example).","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"In order to generate an integration mesh for the interior facets mathcalF_Gamma, we use a new type of Triangulation referred to as SkeletonTriangulation. It can be constructed from a DiscreteModel object as follows:","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"strian = SkeletonTriangulation(model)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"As any other type of Triangulation, an SkeletonTriangulation can be written into a vtk file for its visualization (see next figure, where the interior facets mathcalF_Gamma are clearly observed).","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"writevtk(strian,\"strian\")","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"(Image: )","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"Once we have constructed the triangulations needed in this example, we define the corresponding quadrature rules.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"degree = 2*order\nquad = CellQuadrature(trian,degree)\nbquad = CellQuadrature(btrian,degree)\nsquad = CellQuadrature(strian,degree)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"We still need a way to represent the unit outward normal vector to the boundary partialOmega, and the unit normal vector on the interior faces mathcalF_Gamma. This is done with the get_normal_vector getter.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"nb = get_normal_vector(btrian)\nns = get_normal_vector(strian)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"The get_normal_vector getter takes either a boundary or a skeleton triangulation and returns an object representing the normal vector to the corresponding surface. For boundary triangulations, the returned normal vector is the unit outwards one, whereas for skeleton triangulations the orientation of the returned normal is arbitrary. In the current implementation (Gridap v0.5.0), the unit normal is outwards to the cell with smaller id among the two cells that share an interior facet in mathcalF_Gamma.","category":"page"},{"location":"pages/t006_dg_discretization/#Weak-form-1","page":"6 Poisson equation (with DG)","title":"Weak form","text":"","category":"section"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"With these ingredients we can define the different terms in the weak form. First, we start with the terms a_Omega(cdotcdot) , and b_Omega(cdot) associated with integrals in the volume Omega. This is done as in the tutorial for the Poisson equation.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"a_Ω(u,v) = ∇(v)⊙∇(u)\nb_Ω(v) = v*f\nt_Ω = AffineFETerm(a_Ω,b_Ω,trian,quad)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"The terms a_partialOmega(cdotcdot) and b_partialOmega(cdot) associated with integrals on the boundary partialOmega are defined using an analogous approach. First, we define two functions representing the integrands of the forms a_partialOmega(cdotcdot) and b_partialOmega(cdot). Then, we build an AffineFETerm from these functions and the boundary triangulation and its corresponding quadrature rule:","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"h = L / n\nγ = order*(order+1)\na_∂Ω(u,v) = (γ/h)*v*u - v*(∇(u)⋅nb) - (∇(v)⋅nb)*u\nb_∂Ω(v) = (γ/h)*v*g - (∇(v)⋅nb)*g\nt_∂Ω = AffineFETerm(a_∂Ω,b_∂Ω,btrian,bquad)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"Note that in the definition of the functions a_∂Ω and b_∂Ω, we have used the object nb representing the outward unit normal to the boundary partialOmega. The code definition of a_∂Ω and b_∂Ω is indeed very close to the mathematical definition of the forms a_partialOmega(cdotcdot) and b_partialOmega(cdot).","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"Finally, we need to define the term a_Gamma(cdotcdot) integrated on the interior facets mathcalF_Gamma. In this case, we use a LinearFETerm since the terms integrated on the interior facets only contribute to the system matrix and not to the right-hand-side vector.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"a_Γ(u,v) = (γ/h)*jump(v*ns)⊙jump(u*ns) - jump(v*ns)⊙mean(∇(u)) - mean(∇(v))⊙jump(u*ns)\nt_Γ = LinearFETerm(a_Γ,strian,squad)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"Note that the arguments v, u of function a_Γ represent a test and trial function restricted to the interior facets mathcalF_Gamma. As mentioned before in the presentation of the DG formulation, the restriction of a function vin V to the interior faces leads to two different values v^+ and v^- . In order to compute jumps and averages of the quantities v^+ and v^-, we use the functions jump and mean, which represent the jump and mean value operators lbracklbrack cdot rbrackrbrack and cdot respectively. Note also that we have used the object ns representing the unit normal vector on the interior facets. As a result, the notation used to define function a_Γ is very close to the mathematical definition of the terms in the bilinear form a_Gamma(cdotcdot).","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"Once the different terms of the weak form have been defined, we build and solve the FE problem.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"op = AffineFEOperator(U,V,t_Ω,t_∂Ω,t_Γ)\nuh = solve(op)","category":"page"},{"location":"pages/t006_dg_discretization/#Discretization-error-1","page":"6 Poisson equation (with DG)","title":"Discretization error","text":"","category":"section"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"We end this tutorial by quantifying the discretization error associated with the computed numerical solution uh. In DG methods a simple error indicator is the jump of the computed (discontinuous) approximation on the interior faces. This quantity can be easily computed in Gridap as follows. First, we need to restrict the computed solution uh to the skeleton triangulation.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"uh_Γ = restrict(uh,strian)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"The resulting object uh_Γ is an object which represents the two values u^+_h, u^-_h of the solution u_h restricted to the interior facets mathcalF_Gamma. We compute and visualize the jump of these values as follows (see next figure):","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"writevtk(strian,\"jumps\",cellfields=[\"jump_u\"=>jump(uh_Γ)])","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"Note that the jump of the numerical solution is very small, close to the machine precision (as expected in this example with manufactured solution). (Image: )","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"A more rigorous way of quantifying the error is to measure it with a norm. Here, we use the L^2 and H^1 norms, namely","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"beginaligned\n w _L^2^2 doteq int_Omega w^2 textdOmega \n w _H^1^2 doteq int_Omega w^2 + nabla w cdot nabla w textdOmega\nendaligned","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"The discretization error can be computed in this example as the difference of the manufactured and numerical solutions.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"e = u - uh","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"We compute the error norms as follows. First, we implement the integrands of the norms we want to compute.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"l2(u) = u*u\nh1(u) = a_Ω(u,u) + l2(u)","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"Then, we compute the corresponding integrals with the integrate function.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"el2 = sqrt(sum( integrate(l2(e),trian,quad) ))\neh1 = sqrt(sum( integrate(h1(e),trian,quad) ))","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"The integrate function returns a lazy object representing the contribution to the integral of each cell in the underlying triangulation. To end up with the desired error norms, one has to sum these contributions and take the square root. You can check that the computed error norms are close to machine precision (as one would expect).","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"tol = 1.e-10\n@assert el2 < tol\n@assert eh1 < tol","category":"page"},{"location":"pages/t006_dg_discretization/#References-1","page":"6 Poisson equation (with DG)","title":"References","text":"","category":"section"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"[1] D. N. Arnold, F. Brezzi, B. Cockburn, and L. Donatella Marini. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM Journal on Numerical Analysis, 39 (5):1749–1779, 2001. doi:10.1137/S0036142901384162.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"[2] B. Cockburn, G. Kanschat, and D. Schötzau. An equal-order DG method for the incompressible Navier-Stokes equations. Journal of Scientific Computing, 40(1-3):188–210, 2009. doi:10.1007/s10915-008-9261-1.","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"","category":"page"},{"location":"pages/t006_dg_discretization/#","page":"6 Poisson equation (with DG)","title":"6 Poisson equation (with DG)","text":"This page was generated using Literate.jl.","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"EditURL = \"https://github.com/gridap/Tutorials/blob/master/src/stokes.jl\"","category":"page"},{"location":"pages/t009_stokes/#Tutorial-9:-Stokes-equation-1","page":"9 Stokes equation","title":"Tutorial 9: Stokes equation","text":"","category":"section"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"(Image: ) (Image: )","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"note: Note\nThis tutorial is under construction, but the code below is already functional.","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"Driver that computes the lid-driven cavity benchmark at low Reynolds numbers when using a mixed FE Q(k)/Pdisc(k-1).","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"Load Gridap library","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"using Gridap","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"Discrete model","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"n = 100\ndomain = (0,1,0,1)\npartition = (n,n)\nmodel = CartesianDiscreteModel(domain, partition)","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"Define Dirichlet boundaries","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"labels = get_face_labeling(model)\nadd_tag_from_tags!(labels,\"diri1\",[6,])\nadd_tag_from_tags!(labels,\"diri0\",[1,2,3,4,5,7,8])","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"Define test FESpaces (Q2/P1(disc) pair)","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"V = TestFESpace(\n reffe=:QLagrangian, conformity=:H1, valuetype=VectorValue{2,Float64},\n model=model, labels=labels, order=2, dirichlet_tags=[\"diri0\",\"diri1\"])\nQ = TestFESpace(\n reffe=:PLagrangian, conformity=:L2, valuetype=Float64,\n model=model, order=1, constraint=:zeromean)\nY = MultiFieldFESpace([V,Q])","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"Define trial FESpaces from Dirichlet values","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"u0 = VectorValue(0,0)\nu1 = VectorValue(1,0)\nU = TrialFESpace(V,[u0,u1])\nP = TrialFESpace(Q)\nX = MultiFieldFESpace([U,P])","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"Define integration mesh and quadrature","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"trian = get_triangulation(model); degree = 2\nquad = CellQuadrature(trian,degree)","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"Define and solve the FE problem","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"function a(x,y)\n v,q = y\n u,p = x\n ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u)\nend\nt_Ω = LinearFETerm(a,trian,quad)\nop = AffineFEOperator(X,Y,t_Ω)\nuh, ph = solve(op)","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"Export results to vtk","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"writevtk(trian,\"results\",order=2,cellfields=[\"uh\"=>uh,\"ph\"=>ph])","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"","category":"page"},{"location":"pages/t009_stokes/#","page":"9 Stokes equation","title":"9 Stokes equation","text":"This page was generated using Literate.jl.","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"EditURL = \"https://github.com/gridap/Tutorials/blob/master/src/isotropic_damage.jl\"","category":"page"},{"location":"pages/t010_isotropic_damage/#Tutorial-10:-Isotropic-damage-model-1","page":"10 Isotropic damage model","title":"Tutorial 10: Isotropic damage model","text":"","category":"section"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"(Image: ) (Image: )","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"note: Note\nThis tutorial is under construction, but the code below is already functional.","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"using Gridap\nusing LinearAlgebra","category":"page"},{"location":"pages/t010_isotropic_damage/#Model-definition-1","page":"10 Isotropic damage model","title":"Model definition","text":"","category":"section"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"Elastic branch","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"const E = 3.0e10 # Pa\nconst ν = 0.3 # dim-less\nconst λ = (E*ν)/((1+ν)*(1-2*ν))\nconst μ = E/(2*(1+ν))\n@law σe(ε) = λ*tr(ε)*one(ε) + 2*μ*ε # Pa\nτ(ε) = sqrt(ε ⊙ σe(ε)) # Pa^(1/2)","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"Damage","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"const σ_u = 4.0e5 # Pa\nconst r_0 = σ_u / sqrt(E) # Pa^(1/2)\nconst H = 0.5 # dim-less\n\nfunction d(r)\n 1 - q(r)/r\nend\n\nfunction q(r)\n r_0 + H*(r-r_0)\nend","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"Update of the state variables","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"@law function update(ε_in,r_in,d_in)\n τ_in = τ(ε_in)\n if τ_in <= r_in\n r_out = r_in\n d_out = d_in\n damaged = false\n else\n r_out = τ_in\n d_out = d(r_out)\n damaged = true\n end\n damaged, r_out, d_out\nend","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"Constitutive law and its linearization","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"@law function σ(ε_in,r_in,d_in)\n\n _, _, d_out = update(ε_in,r_in,d_in)\n\n (1-d_out)*σe(ε_in)\n\nend\n\n@law function dσ(dε_in,ε_in,state)\n\n damaged, r_out, d_out = state\n\n if ! damaged\n return (1-d_out)*σe(dε_in)\n\n else\n c_inc = ((q(r_out) - H*r_out)*(σe(ε_in) ⊙ dε_in))/(r_out^3)\n return (1-d_out)*σe(dε_in) - c_inc*σe(ε_in)\n\n end\nend","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"max dead load","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"const b_max = VectorValue(0.0,0.0,-(9.81*2.5e3))","category":"page"},{"location":"pages/t010_isotropic_damage/#L2-projection-1","page":"10 Isotropic damage model","title":"L2 projection","text":"","category":"section"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"form Gauss points to a Lagrangian piece-wise discontinuous space","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"function project(q,trian,quad,order)\n\n a(u,v) = u*v\n l(v) = v*q\n t_Ω = AffineFETerm(a,l,trian,quad)\n\n V = TestFESpace(\n reffe=:Lagrangian, valuetype=Float64, order=order,\n triangulation=trian, conformity=:L2)\n\n U = TrialFESpace(V)\n op = AffineFEOperator(U,V,t_Ω)\n qh = solve(op)\n qh\nend","category":"page"},{"location":"pages/t010_isotropic_damage/#Main-function-1","page":"10 Isotropic damage model","title":"Main function","text":"","category":"section"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"function main(;n,nsteps)\n\n r = 12\n domain = (0,r,0,1,0,1)\n partition = (r*n,n,n)\n model = CartesianDiscreteModel(domain,partition)\n\n labeling = get_face_labeling(model)\n add_tag_from_tags!(labeling,\"supportA\",[1,3,5,7,13,15,17,19,25])\n add_tag_from_tags!(labeling,\"supportB\",[2,4,6,8,14,16,18,20,26])\n add_tag_from_tags!(labeling,\"supports\",[\"supportA\",\"supportB\"])\n\n order = 1\n\n V = TestFESpace(\n reffe=:Lagrangian, valuetype=VectorValue{3,Float64}, order=order,\n model=model,\n labels = labeling,\n dirichlet_tags=[\"supports\"])\n\n U = TrialFESpace(V)\n\n trian = Triangulation(model)\n degree = 2*order\n quad = CellQuadrature(trian,degree)\n\n r = QPointCellField(r_0,trian,quad)\n d = QPointCellField(0.0,trian,quad)\n\n nls = NLSolver(show_trace=true, method=:newton)\n solver = FESolver(nls)\n\n function step(uh_in,factor)\n\n b = factor*b_max\n res(u,v) = ( ε(v) ⊙ σ(ε(u),r,d) ) - v⋅b\n function jac(u,du,v)\n state = update(ε(u),r,d)\n ε(v) ⊙ dσ(ε(du),ε(u),state)\n end\n t_Ω = FETerm(res,jac,trian,quad)\n op = FEOperator(U,V,t_Ω)\n\n uh_out, = solve!(uh_in,solver,op)\n\n update_state_variables!(quad,update,ε(uh_out),r,d)\n\n uh_out\n end\n\n factors = collect(1:nsteps)*(1/nsteps)\n uh = zero(V)\n\n for (istep,factor) in enumerate(factors)\n\n println(\"\\n+++ Solving for load factor $factor in step $istep of $nsteps +++\\n\")\n\n uh = step(uh,factor)\n dh = project(d,trian,quad,order)\n rh = project(r,trian,quad,order)\n\n writevtk(\n trian,\"results_$(lpad(istep,3,'0'))\",\n cellfields=[\"uh\"=>uh,\"epsi\"=>ε(uh),\"damage\"=>dh,\n \"threshold\"=>rh,\"sigma_elast\"=>σe(ε(uh))])\n\n end\n\nend","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"Run!","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"main(n=6,nsteps=20)","category":"page"},{"location":"pages/t010_isotropic_damage/#Results-1","page":"10 Isotropic damage model","title":"Results","text":"","category":"section"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"Animation of the load history using for main(n=8,nsteps=30) (Image: )","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"","category":"page"},{"location":"pages/t010_isotropic_damage/#","page":"10 Isotropic damage model","title":"10 Isotropic damage model","text":"This page was generated using Literate.jl.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"EditURL = \"https://github.com/gridap/Tutorials/blob/master/src/fsi_tutorial.jl\"","category":"page"},{"location":"pages/t011_fsi_tutorial/#Tutorial-11:-Fluid-Structure-Interaction-1","page":"11 Fluid-Structure Interaction","title":"Tutorial 11: Fluid-Structure Interaction","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"(Image: ) (Image: )","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"In this tutorial, we will learn","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"How to solve a surface-coupled multi-physics problem.\nConstruct FE spaces defined in different domains.\nDefine interface triangulations and integrate on them.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Problem Statement\nStrong form\nGeometry and Discrete model\nBoundary conditions and properties\nNumerical scheme\nFE spaces\nWeak form\nNumerical integration\nAlgebraic system of equations\nPost-processing\nVisualization\nQuantities of interest","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"probStat\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#Problem-statement-1","page":"11 Fluid-Structure Interaction","title":"Problem statement","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"strongForm\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#Strong-form-1","page":"11 Fluid-Structure Interaction","title":"Strong form","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"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 FD and Gamma_rm FN the fluid boundaries with Dirichlet and Neumann conditions, respectively. The Fluid-Structure Interaction (FSI) problem reads:","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"find u_rm F, p_rm F and u_rm S such that","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"leftlbrace\nbeginaligned\n-nablacdotboldsymbolsigma_rm F = f text in Omega_rm F\nnablacdot u_rm F = 0 text in Omega_rm F\n-nablacdotboldsymbolsigma_rm S = s text in Omega_rm S\nendaligned\nright","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"satisfying the Dirichlet and Neumann boundary conditions","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"leftlbrace\nbeginaligned\nu_rm F = g_rm F text on Gamma_rm FD\nn_rm F cdot boldsymbolsigma_rm F = 0 text on Gamma_rm FN\nu_rm S = g_rm S text on Gamma_rm SD\nendaligned\nright","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"and the kinematic and dynamic conditions at the fluid-solid interface","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"leftlbrace\nbeginaligned\nu_rm F = u_rm S text on Gamma_rm FS\nn_rm F cdot boldsymbolsigma_rm F + n_rm S cdot boldsymbolsigma_rm S = 0 text on Gamma_rm FS\nendaligned\nright","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Where boldsymbolsigma_rm F(u_rm Fp_rm F)=2mu_rm Fboldsymbolvarepsilon(u_rm F) - p_rm FmathbfI and boldsymbolsigma_rm S(u_rm S)=2mu_rm Sboldsymbolvarepsilon(u_rm S) +lambda_rm Str(boldsymbolvarepsilon(u_rm S))mathbfI.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"geometry\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#Geometry-and-Discrete-model-1","page":"11 Fluid-Structure Interaction","title":"Geometry and Discrete model","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"In this tutorial we solve the benchmark descrived in [1], consisting on a flow over an elastic flag after a cylinder. The computational domain is defined by a channel of size Omega doteq (045)times(0041), with an embedded cylinder of radius R=005 and center at C=(0202). The associated FE triangulation is denoted by mathcalT, the fluid and solid domain and their associated triangulations will be denoted by Omega_rm F, Omega_rm S, mathcalT_rm F and mathcalT_rm S, respectively.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"In order to load the discrete model we first setup Gridap","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"using Gridap","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"The discrete model for the elastic flag problem is generated by loading the \"../models/elasticFlag.json\" file.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"model = DiscreteModelFromFile(\"../models/elasticFlag.json\")","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"We can inspect the loaded geometry and associated parts by printing to a vtk file:","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"writevtk(model,\"model\")","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"This will produce an output in which we can identify the different parts of the domain, with the associated labels and tags.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Part Notation Label Tag\nSolid-cylinder wall Gamma_rm SD_cyl \"fixed\" 1\nFluid-solid interface Gamma_rm FS \"interface\" 2\nChannel inlet Gamma_rm FD_in \"inlet\" 3\nChannel outlet Gamma_rm FN_out \"outlet\" 4\nChannel walls Gamma_rm FD_wall \"noSlip\" 5\nFluid-cylinder wall Gamma_rm FD_cyl \"cylinder\" 6\nFluid domain Omega_rm F \"fluid\" 7\nSolid domain Omega_rm S \"solid\" 8","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"(Image: )","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"conditions\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#External-conditions-and-properties-1","page":"11 Fluid-Structure Interaction","title":"External conditions and properties","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#Boundary-conditions-1","page":"11 Fluid-Structure Interaction","title":"Boundary conditions","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"We apply Dirichlet boundary conditions at the channel inlet, upper and lower boundaries and on the cylinder. A parabolic profile is enforced at the channel inlet, while a no-slip condition is imposed on the remaining Dirichlet boundaries.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"leftlbrace\nbeginaligned\nu_rm Fin(xy) = 15Ufracy(H y)left(fracH2right)^2quadtextrmon Gamma_rm FD_in\nu_rm F0(xy) = 0 0quadtextrmon Gamma_rm FD_wallcupGamma_rm FD_cyl\nu_rm S0(xy) = 0 0quadtextrmon Gamma_rm SD_cyl\nendaligned\nright","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"const U = 1.0\nconst H = 0.41\nuf_in(x) = VectorValue( 1.5 * U * x[2] * ( H - x[2] ) / ( (H/2)^2 ), 0.0 )\nuf_0(x) = VectorValue( 0.0, 0.0 )\nus_0(x) = VectorValue( 0.0, 0.0 )","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"We consider a free tranction condition at the channel outlet","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"n_rm F cdot boldsymbolsigma_rm F = mathbf0quadtextrmon Gamma_rm FN","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"hN(x) = VectorValue( 0.0, 0.0 )\np_jump(x) = 0.0","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"forces\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#External-forces-1","page":"11 Fluid-Structure Interaction","title":"External forces","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"In this test, the body forces acting on the fluid an solid are zero.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"f(x) = VectorValue( 0.0, 0.0 )\ns(x) = VectorValue( 0.0, 0.0 )\ng(x) = 0.0","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"properties\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#Material-properties-1","page":"11 Fluid-Structure Interaction","title":"Material properties","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"We use a linear elastic constitutive law for the elastic flag. Given the Young's modulus E and the Poisson ratio nu, we can compute the Lamé constants, lambda and mu, using the following function:","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"function lame_parameters(E,ν)\n λ = (E*ν)/((1+ν)*(1-2*ν))\n μ = E/(2*(1+ν))\n (λ, μ)\nend","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Then, we get the Lamé parameters for a solid with E=10 MPa and nu=033.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"const E_s = 1.0\nconst ν_s = 0.33\nconst (λ_s,μ_s) = lame_parameters(E_s,ν_s)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"The Cauchy stress tensor for the solid part is defined by sigma_s = 2muvarepsilon + lambda tr(varepsilon)mathbfI. Note that we use the trace operator from the LinearAlgebra package. With the macro @law we are able to define a function whose arguments depend on the coordinates, without the need of passing such coordinates as an argument.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"using LinearAlgebra: tr\n@law σ_s(ε) = λ_s*tr(ε)*one(ε) + 2*μ_s*ε","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"For the fluid part, we only need to define the viscosity mu_f, which we set to 001.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"const μ_f = 0.01","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"The Cauchy stress tensor for the fluid part is given by sigma_f = sigma^dev_f - pmathbfI, with sigma^dev_f=2mu_f the deviatoric part of the stress. Since we use a mixed form with the pressure p as an unknown, the stress law will only involve the deviatoric part.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"@law σ_dev_f(ε) = 2*μ_f*ε","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"numericalScheme\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#Numerical-scheme-1","page":"11 Fluid-Structure Interaction","title":"Numerical scheme","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"feSpace\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#FE-Spaces-1","page":"11 Fluid-Structure Interaction","title":"FE Spaces","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"In this tutorial we use an inf-sup stable velocity-pressure pair, P_kP_k-1 elements, with continuous velocities and pressures. We select the same velocity interpolation space for the fluid and the solid parts, defined as follows","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"V_rm X doteq v in H^1(Omega_rm X)^d v_Tin P_k(T)^d text for all TinmathcalT_rm X ","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"where T denotes an arbitrary cell of the FE mesh mathcalT_rm X, and P_k(T) is the usual Lagrangian FE space of order k defined on a mesh of triangles or tetrahedra. Note that the sub-index (cdot)_rm X stands for the fluid or solid parts, (cdot)_rm F or (cdot)_rm S, respectively. On the other hand, the space for the pressure is only defined in the fluid domain, Omega_rm F, and is given by","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Q doteq q in C^0(Omega) q_Tin P_k-1(T) text for all TinmathcalT_rm F","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Before creating the FE spaces, we first need to define the discrete model of each of the sub-domains, which are constructed restricting the global discrete model to the corresponding part. This is done by calling the DiscreteModel function with the desired geometrical part label, i.e. \"solid\" and \"fluid\", respectively.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"model_solid = DiscreteModel(model,\"solid\")\nmodel_fluid = DiscreteModel(model,\"fluid\")","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"In what follows we will assume a second-order veloticty interpolation, i.e. k=2","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"k = 2","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Having set up all the ingredients, we can create the different FE spaces for the test functions. For the velocity FE spaces we call the TestFESpace function with the corresponding discrete model, using a 2-dimensional VectorValue type, a :Lagrangian reference FE element of order k and conformity :H1. Note that we assign different Dirichlet boundary labels for the two different parts, generating the variational spaces with homogeneous Dirichlet boundary conditions, V_rm F0 and V_rm S0 .","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Vf = TestFESpace(\n model=model_fluid,\n valuetype=VectorValue{2,Float64},\n reffe=:Lagrangian,\n order=k,\n conformity =:H1,\n dirichlet_tags=[\"inlet\", \"noSlip\", \"cylinder\"])\n\nVs = TestFESpace(\n model=model_solid,\n valuetype=VectorValue{2,Float64},\n reffe=:Lagrangian,\n order=k,\n conformity =:H1,\n dirichlet_tags=[\"fixed\"])","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"For the pressure test FE space, we use the fluid discrete model, a scalar value type, a :Lagrangian reference FE element of order k-1 and :C0 conformity.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Qf = TestFESpace(\n model=model_fluid,\n valuetype=Float64,\n order=k-1,\n reffe=:Lagrangian,\n conformity=:C0)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"The trial FE spaces are generated from the test FE spaces, adding the corresponding function for the various Dirichlet boundaries, leading to U_rm Fg_rm F, U_rm Sg_rm S and P_rm F.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Uf = TrialFESpace(Vf,[uf_in, uf_0, uf_0])\nUs = TrialFESpace(Vs,[us_0])\nPf = TrialFESpace(Qf)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Finally, we glue the test and trial FE spaces together, defining a unique test and trial space for all the fields using the MultiFieldFESpace function. That is Y=V_rm S0 V_rm F0 Q_rm F^T and X=U_rm Sg_rm S U_rm Fg_rm F P_rm F^T","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Y = MultiFieldFESpace([Vs,Vf,Qf])\nX = MultiFieldFESpace([Us,Uf,Pf])","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"weakForm\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#Weak-form-1","page":"11 Fluid-Structure Interaction","title":"Weak form","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"We now introduce the solution and test function vectors as mathbfx^h=mathbfu^h_rm Smathbfu^h_rm F p^h_rm F^T and mathbfy^h=mathbfv^h_rm Smathbfv^h_rm F q^h_rm F^T. The weak form of the coupled FSI problem using the Nitche's method, see for instance [2], reads: find mathbfx^h in X such that","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"a_s(mathbfx^hmathbfy^h) + a_f(mathbfx^hmathbfy^h) + a_fs(mathbfx^hmathbfy^h) = l_s(mathbfy^h) + l_f(mathbfy^h) + l_fGamma_N(mathbfy^h)qquadforallmathbfy^hin Y","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"with the following definitions:","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"a_s(mathbfx^hmathbfy^h)\nis the bilinear form associated with the solid counterpart, defined as","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"a_s(mathbfx^hmathbfy^h)=(varepsilon(mathbfv^h_rm S)boldsymbolsigma_s(mathbfu^h_rm S))","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"function a_s(x,y)\n us,uf,p = x\n vs,vf,q = y\n\tε(vs) ⊙ σ_s(ε(us))\nend","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"a_f(mathbfx^hmathbfy^h)\nis the bilinear form associated with the fluid counterpart, defined as","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"a_f(mathbfx^hmathbfy^h)=(varepsilon(mathbfv^h_rm F)boldsymbolsigma^dev_f(mathbfu^h_rm F)) - (nablacdotmathbfv^h_rm Fp_rm F) + (q_rm F nablacdotmathbfu^h_rm F)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"function a_f(x,y)\n us,uf,p = x\n vs,vf,q = y\n (ε(vf) ⊙ σ_dev_f(ε(uf))) - (∇⋅vf)*p + q*(∇⋅uf)\nend","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"a_fs(mathbfx^hmathbfy^h)\nis the bilinear form associated with the coupling between fluid and solid counterparts. To difine this form we use the well known Nitsche's method, which enforces the continuity of fluid and solid velocities as well as the continuity of the normal stresses, see for instance [2]. The final expression for this term reads:","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"beginaligned\na_fs(mathbfx^hmathbfy^h)=langlegammafracmu_fh(mathbfv^h_rm F-mathbfv^h_rm S)(mathbfu^h_rm F-mathbfu^h_rm S)rangle_Gamma_rm FS\n- langle(mathbfv^h_rm F-mathbfv^h_rm S)mathbfn_rm FScdotboldsymbolsigma^dev_f(mathbfu^h_rm F)rangle_Gamma_rm FS\n+ langle(mathbfv^h_rm F-mathbfv^h_rm S)p^h_rm Fmathbfn_rm FSrangle_Gamma_rm FS\n- chilanglemathbfn_rm FScdotboldsymbolsigma^dev_f(mathbfv^h_rm F)(mathbfu^h_rm F-mathbfu^h_rm S)rangle_Gamma_rm FS\n+ langle q^h_rm Fmathbfn_rm FS(mathbfu^h_rm F-mathbfu^h_rm S)rangle_Gamma_rm FS\nendaligned","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Where chi is a parameter that can take values 10 or -10 and it is used to define the symmetric or antisymmetric version of the method, respectively.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"const γ = 1.0\nconst χ = -1.0\nfunction a_fs(x,y)\n us_Γ, uf_Γ, p_Γ = x\n vs_Γ, vf_Γ, q_Γ = y\n uf = jump(uf_Γ)\n p = jump(p_Γ)\n us = -jump(us_Γ)\n vf = jump(vf_Γ)\n q = jump(q_Γ)\n vs = -jump(vs_Γ)\n εuf = jump(ε(uf_Γ))\n εvf = jump(ε(vf_Γ))\n penaltyTerms = α*vf⋅uf - α*vf⋅us - α*vs⋅uf + α*vs⋅us\n\tintegrationByParts = ( vf⋅(p*n_Γfs) - vf⋅(n_Γfs⋅σ_dev_f(εuf)) ) - ( vs⋅(p*n_Γfs) - vs⋅(n_Γfs⋅σ_dev_f(εuf)) )\n\tsymmetricTerms = ( q*(n_Γfs⋅uf) - χ*(n_Γfs⋅σ_dev_f(εvf))⋅uf ) - ( q*(n_Γfs⋅us) - χ*(n_Γfs⋅σ_dev_f(εvf))⋅us )\n penaltyTerms + integrationByParts + symmetricTerms\nend","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"l_s(mathbfy^h)\nis the linear form associated with the solid counterpart, defined as","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"l_s(mathbfy^h)=(mathbfv^h_rm Ss)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"function l_s(y)\n vs,vf,q = y\n vs⋅s\nend","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"l_f(mathbfy^h)\nis the linear form associated with the fluid counterpart, defined as","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"l_f(mathbfy^h)=(mathbfv^h_rm Ff)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"function l_f(y)\n vs,vf,q = y\n vf⋅f + q*g\nend","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"l_fGamma_N(mathbfy^h)\nis the linear form associated with the fluid Neumann boundary condition, defined as","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"l_fGamma_N(mathbfy^h)=langlemathbfv^h_rm Fhrangle_Gamma_N","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"function l_f_Γn(y)\n vs,vf,q = y\n vf⋅hN\nend","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"integration\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#Numerical-integration-1","page":"11 Fluid-Structure Interaction","title":"Numerical integration","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"To define the quadrature rules used in the numerical integration of the different terms, we first need to generate the domain triangulation. Here we create the triangulation of the global domain, mathcalT.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"trian = Triangulation(model)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"The solid and fluid triangulations, mathcalT_rm F and mathcalT_rm S, are constructed from the discrete models restricted to the respective parts.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"trian_solid = Triangulation(model_solid)\ntrian_fluid = Triangulation(model_fluid)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"We also generate the triangulation and associated outer normal field for the outlet boundary, Gamma_rm FN_out, which will be used to impose a Neumann condition.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"trian_Γout = BoundaryTriangulation(model,\"outlet\")\nn_Γout = get_normal_vector(trian_Γout)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Finally, to impose the interface condition between solid and fluid, we will need the triangulation and normal field of such interface, Gamma_rm FS. The interface triangulation is generated by calling the InterfaceTriangulation function specifying the two interfacing domain models. Note that the normal field will point outwards with respect to the first entry.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"trian_Γfs = InterfaceTriangulation(model_fluid,model_solid)\nn_Γfs = get_normal_vector(trian_Γfs)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"From the interface triangulation we can obtain the interface elements length, h, and the penalty parameter, alpha=gammafracmu_fh, used in the Nitsche's terms.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"using Gridap.Arrays\nusing LinearAlgebra: norm\nxe = get_cell_coordinates(trian_Γfs)\nhe = apply(x->norm(x[2]-x[1]),xe)\nα = apply(h->γ*μ_f/h, he)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Once we have all the triangulations, we can generate the quadrature rules to be applied each domain.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"degree = 2*k\nquad_solid = CellQuadrature(trian_solid,degree)\nquad_fluid = CellQuadrature(trian_fluid,degree)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Idem for the boundary part.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"bdegree = 2*k\nquad_Γout = CellQuadrature(trian_Γout,bdegree)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Idem for the interface part.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"idegree = 2*k\nquad_Γfs = CellQuadrature(trian_Γfs,idegree)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"algebraic\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#Algebraic-System-of-Equations-1","page":"11 Fluid-Structure Interaction","title":"Algebraic System of Equations","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"After defining the weak form of the problem and the integration quadrature rules to perform the numerical integration, we are ready to assemble the linear system of equations. In this case, the system will have the following structure:","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"beginbmatrix\nmathbfA_urm SmathbfA_urm SF0\nmathbfA_urm FSmathbfA_urm FmathbfA_uprm F\n0mathbfA_uprm FmathbfA_prm F\nendbmatrixbeginbmatrix\nmathbfU^h_rm S\nmathbfU^h_rm F\nmathbfP^h_rm F\nendbmatrix = beginbmatrix\nmathbfF^h_urm S\nmathbfF^h_urm F\nmathbfF^h_prm F\nendbmatrix","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"In order to construct the system we first define the diferent discrete terms using the functions AffineFETerm, that assemble contributions on the left-hand side and the right-hand side of the system, for the fluid and solid interior terms evaluated in the respective triangulations, mathcalT_rm F and mathcalT_rm S. The fluid Neumann boundary condition is assembled using the function FESource, which only affects the right-hand side of the system, evaluating the corresponding form on the triangulation of the outlet boundary, Gamma_rm FN_out. The coupling terms are evaluated at the interface Gamma_rm FS using the function LinearFETerm, assembling only to the left-hand side of the system.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"t_Ω_s= AffineFETerm(a_s,l_s,trian_solid,quad_solid)\nt_Ω_f = AffineFETerm(a_f,l_f,trian_fluid,quad_fluid)\nt_Γfs = LinearFETerm(a_fs,trian_Γfs,quad_Γfs)\nt_Γn_f = FESource(l_f_Γn,trian_Γout,quad_Γout)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"The final FE operator is constructed using the function AffineFEOperator and passing as arguments the trial and test FE spaces, X and Y, and all the different terms previously defined.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"op = AffineFEOperator(X,Y,t_Ω_s,t_Ω_f,t_Γn_f,t_Γfs)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"Finally, we call solve to obtain the solution vector of nodal values mathbfU^h_rm SmathbfU^h_rm FmathbfP^h_rm F^T","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"uhs, uhf, ph = solve(op)","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"postprocess\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#Post-processing-1","page":"11 Fluid-Structure Interaction","title":"Post-processing","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"visualization\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#Visualization-1","page":"11 Fluid-Structure Interaction","title":"Visualization","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"The solution fields mathbfU^h_rm SmathbfU^h_rm FmathbfP^h_rm F^T are defined over all the domain, extended with zeros on the inactive part. Calling the function writevtk passing the global triangulation, we will output the global fields.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"writevtk(trian,\"trian\", cellfields=[\"uhs\" => uhs, \"uhf\" => uhf, \"ph\" => ph])","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"(Image: )","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"However, we can also restrict the fields to the active part by calling the function restrict with the field along with the respective active triangulation.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"uhs_solid = restrict(uhs, trian_solid)\nuhf_fluid = restrict(uhf, trian_fluid)\nph_fluid = restrict(ph, trian_fluid)\nwritevtk(trian_solid,\"trian_solid\",cellfields=[\"uhs\"=>uhs_solid])\nwritevtk(trian_fluid,\"trian_fluid\",cellfields=[\"ph\"=>ph_fluid,\"uhf\"=>uhf_fluid])","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"(Image: )","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"<a name=\"QOIs\"></a>","category":"page"},{"location":"pages/t011_fsi_tutorial/#Quantities-of-Interest-1","page":"11 Fluid-Structure Interaction","title":"Quantities of Interest","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"trian_ΓS = BoundaryTriangulation(model,[\"cylinder\",\"interface\"])\nquad_ΓS = CellQuadrature(trian_ΓS,bdegree)\nn_ΓS = get_normal_vector(trian_ΓS)\nuh_ΓS = restrict(uhf_fluid,trian_ΓS)\nph_ΓS = restrict(ph_fluid,trian_ΓS)\nFD, FL = sum( integrate( (n_ΓS⋅σ_dev_f(ε(uh_ΓS))) - ph_ΓS*n_ΓS, trian_ΓS, quad_ΓS ) )\nprintln(\"Drag force: \", FD)\nprintln(\"Lift force: \", FL)","category":"page"},{"location":"pages/t011_fsi_tutorial/#References-1","page":"11 Fluid-Structure Interaction","title":"References","text":"","category":"section"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"[1] Turek, S., Hron, J., Madlik, M., Razzaq, M., Wobker, H., & Acker, J. F. (2011).* Numerical simulation and benchmarking of a monolithic multigrid solver for fluid-structure interaction problems with application to hemodynamics. In Fluid Structure Interaction II (pp. 193-220). Springer, Berlin, Heidelberg.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"[2] Burman, E., and Fernández, M. A. Stabilized explicit coupling for fluid–structure interaction using Nitsche's method. Comptes Rendus Mathematique 345.8 (2007): 467-472.","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"","category":"page"},{"location":"pages/t011_fsi_tutorial/#","page":"11 Fluid-Structure Interaction","title":"11 Fluid-Structure Interaction","text":"This page was generated using Literate.jl.","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"EditURL = \"https://github.com/gridap/Tutorials/blob/master/src/p_laplacian.jl\"","category":"page"},{"location":"pages/t004_p_laplacian/#Tutorial-4:-p-Laplacian-1","page":"4 p-Laplacian","title":"Tutorial 4: p-Laplacian","text":"","category":"section"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"(Image: ) (Image: )","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"In this tutorial, we will learn","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"How to solve a simple nonlinear PDE in Gridap\nHow to define the weak residual and its Jacobian\nHow to setup and use a nonlinear solver\nHow to define new boundaries from a given discrete model","category":"page"},{"location":"pages/t004_p_laplacian/#Problem-statement-1","page":"4 p-Laplacian","title":"Problem statement","text":"","category":"section"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"The goal of this tutorial is to solve a nonlinear PDE in Gridap. For the sake of simplicity, we consider the p-Laplacian equation as the model problem. Specifically, the PDE we want to solve is: find the scalar-field u such that","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"leftlbrace\nbeginaligned\n-nabla cdot left( nabla u^p-2 nabla u right) = f textin Omega\nu = 0 texton Gamma_0\nu = g texton Gamma_g\nleft( nabla u^p-2 nabla u right)cdot n = 0 texton Gamma_rm N\nendaligned\nright","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"with p2. The computational domain Omega is the one depicted in next figure, which is the same as in the first tutorial. However, we slightly change the boundary conditions here. We impose homogeneous Dirichlet and homogeneous Neumann boundary conditions on Gamma_0 and Gamma_rm N respectively, and in-homogeneous Dirichlet conditions on Gamma_g. The Dirichlet boundaries Gamma_0 and Gamma_g are defined as the closure of the green and blue surfaces in next figure respectively, whereas the Neumann boundary is Gamma_rm NdoteqpartialOmega setminus (Gamma_0cupGamma_g). In this example, we consider the values p=3, f=1, and g=2.","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"(Image: )","category":"page"},{"location":"pages/t004_p_laplacian/#Numerical-scheme-1","page":"4 p-Laplacian","title":"Numerical scheme","text":"","category":"section"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"We discretize the problem with conforming Lagrangian FE spaces. For this formulation, the nonlinear weak form reads: find uin U_g such that r(u)(v) = 0 for all vin V_0. As in previous tutorials, the space U_g is the set of functions in H^1(Omega) that fulfill the Dirichlet boundary conditions, whereas V_0 is composed by functions in H^1(Omega) that vanish at the Dirichlet boundary. The weak residual r(u) evaluated at a function uin U_g is the linear form defined as","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"r(u)(v) doteq int_Omega nabla v cdot left( nabla u^p-2 nabla u right) rm dOmega - int_Omega v f rm dOmega","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"In order to solve this nonlinear weak equation, we consider a Newton-Raphson method, which is associated with a linearization of the problem in an arbitrary direction delta uin V_0, namely r(u+delta u)(v)approx r(u)(v) + j(u)(delta uv). In previous formula, j(u) is the Jacobian evaluated at uin U_g, which is the bilinear form","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"j(u)(delta uv) = int_Omega nabla v cdot left( nabla u^p-2 nabla delta u right) rm dOmega + (p-2) int_Omega nabla v cdot left( nabla u^p-4 (nabla u cdot nabla delta u) nabla u right) rm dOmega","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"Note that the solution of this nonlinear PDE with a Newton-Raphson method, will require to discretize both the residual r and the Jacobian j. In Gridap, this is done by following an approach similar to the one already shown in previous tutorials for discretizing the bilinear and linear forms associated with a linear FE problem. The specific details are discussed now.","category":"page"},{"location":"pages/t004_p_laplacian/#Discrete-model-1","page":"4 p-Laplacian","title":"Discrete model","text":"","category":"section"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"As in previous tutorials, the first step to solve the PDE is to load a discretization of the computational domain. In this case, we load the model from the same file as in the first tutorial","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"using Gridap\nmodel = DiscreteModelFromFile(\"../models/model.json\")","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"As stated before, we want to impose Dirichlet boundary conditions on Gamma_0 and Gamma_g, but none of these boundaries is identified in the model. E.g., you can easily see by writing the model in vtk format","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"writevtk(model,\"model\")","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"and by opening the file \"model_0\" in Paraview that the boundary identified as \"sides\" only includes the vertices in the interior of Gamma_0, but here we want to impose Dirichlet boundary conditions in the closure of Gamma_0, i.e., also on the vertices on the contour of Gamma_0. Fortunately, the objects on the contour of Gamma_0 are identified with the tag \"sides_c\" (see next figure). Thus, the Dirichlet boundary Gamma_0 can be built as the union of the objects identified as \"sides\" and \"sides_c\".","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"(Image: )","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"Gridap provides a convenient way to create new object identifiers (referred to as \"tags\") from existing ones. First, we need to extract from the model, the object that holds the information about the boundary identifiers (referred to as FaceLabeling):","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"labels = get_face_labeling(model)","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"Then, we can add new identifiers (aka \"tags\") to it. In the next line, we create a new tag called \"diri0\" as the union of the objects identified as \"sides\" and \"sides_c\", which is precisely what we need to represent the closure of the Dirichlet boundary Gamma_0.","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"add_tag_from_tags!(labels,\"diri0\",[\"sides\", \"sides_c\"])","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"We follow the same approach to build a new identifier for the closure of the Dirichlet boundary Gamma_g. In this case, the boundary is expressed as the union of the objects identified with the tags \"circle\", \"circle_c\", \"triangle\", \"triangle_c\", \"square\", \"square_c\". Thus, we create a new tag for Gamma_g, called \"dirig\" simply as follows:","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"add_tag_from_tags!(labels,\"dirig\",\n [\"circle\",\"circle_c\", \"triangle\", \"triangle_c\", \"square\", \"square_c\"])","category":"page"},{"location":"pages/t004_p_laplacian/#FE-Space-1","page":"4 p-Laplacian","title":"FE Space","text":"","category":"section"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"Now, we can build the FE space by using the newly defined boundary tags.","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"V0 = TestFESpace(\n reffe=:Lagrangian, order=1, valuetype=Float64,\n conformity=:H1, model=model, labels=labels,\n dirichlet_tags=[\"diri0\", \"dirig\"])","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"The construction of this space is essentially the same as in the first tutorial (we build a continuous scalar-valued Lagrangian interpolation of first order). However, we also pass here the labels object (that contains the newly created boundary tags). From this FE space, we define the trial FE spaces","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"g = 1\nUg = TrialFESpace(V0,[0,g])","category":"page"},{"location":"pages/t004_p_laplacian/#Nonlinear-FE-problem-1","page":"4 p-Laplacian","title":"Nonlinear FE problem","text":"","category":"section"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"At this point, we are ready to build the nonlinear FE problem. To this end, we need to define the weak residual and also its corresponding Jacobian. This is done following a similar procedure to the one considered in previous tutorials to define the bilinear and linear forms associated with linear FE problems. In this case, instead of an AffineFETerm (which is for linear problems), we use a FETerm which accounts for the non-linear case. An instance of FETerm is constructed by providing the integrands of the weak residual and its Jacobian (in a similar way an AffineFETerm is constructed from the integrands of the bilinear and linear forms).","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"On the one hand, the integrand of the weak residual is built as follows","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"using LinearAlgebra: norm\nconst p = 3\n@law flux(∇u) = norm(∇u)^(p-2) * ∇u\nf(x) = 1\nres(u,v) = ∇(v)⊙flux(∇(u)) - v*f","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"Function res is the one representing the integrand of the weak residual r(u)(v). The first argument of function res stands for the function uin U_g, where the residual is evaluated, and the second argument stands for a generic test function vin V_0. Note that we have used the macro @law to construct the \"constitutive law\" that relates the nonlinear flux with the gradient of the solution.","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"On the other hand, we implement a function jac representing the integrand of the Jacobian","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"@law dflux(∇du,∇u) =\n (p-2)*norm(∇u)^(p-4)*(∇u ⊙ ∇du)*∇u + norm(∇u)^(p-2) * ∇du\njac(u,du,v) = ∇(v)⊙dflux(∇(du),∇(u))","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"The first argument of function jac stands for function uin U_g, where the Jacobian is evaluated. The second argument is a test function vin V_0, and the third argument represents an arbitrary direction delta u in V_0. Note that we have also used the macro @law to define the linearization of the nonlinear flux.","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"With these functions, we build the FETerm as follows:","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"trian = Triangulation(model)\ndegree=2\nquad = CellQuadrature(trian,degree)\nt_Ω = FETerm(res,jac,trian,quad)","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"We build the FETerm by passing in the first and second arguments the functions that represent the integrands of the residual and Jacobian respectively. The other two arguments, are the triangulation and quadrature used to perform the integrals numerically. From this FETerm object, we finally construct the nonlinear FE problem","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"op = FEOperator(Ug,V0,t_Ω)","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"Here, we have constructed an instance of FEOperator, which is the type that represents a general nonlinear FE problem in Gridap. The constructor takes the test and trial spaces, and the FETerms objects describing the corresponding weak form (in this case only a single term).","category":"page"},{"location":"pages/t004_p_laplacian/#Nonlinear-solver-phase-1","page":"4 p-Laplacian","title":"Nonlinear solver phase","text":"","category":"section"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"We have already built the nonlinear FE problem. Now, the remaining step is to solve it. In Gridap, nonlinear (and also linear) FE problems can be solved with instances of the type FESolver.","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"We construct an instance of FESolver as follows:","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"using LineSearches: BackTracking\nnls = NLSolver(\n show_trace=true, method=:newton, linesearch=BackTracking())\nsolver = FESolver(nls)","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"Note that the NLSolver function used above internally calls the nlsolve function of the NLsolve package with the provided key-word arguments. Thus, one can use any of the nonlinear methods available via the function nlsolve to solve the nonlinear FE problem. Here, we have selected a Newton-Raphson method with a back-tracking line-search from the LineSearches package.","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"We are finally in place to solve the nonlinear FE problem. The initial guess is a FEFunction, which we build from a vector of random (free) nodal values:","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"import Random\nRandom.seed!(1234)\nx = rand(Float64,num_free_dofs(Ug))\nuh0 = FEFunction(Ug,x)\nuh, = solve!(uh0,solver,op)","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"We finish this tutorial by writing the computed solution for visualization (see next figure).","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"writevtk(trian,\"results\",cellfields=[\"uh\"=>uh])","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"(Image: )","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"","category":"page"},{"location":"pages/t004_p_laplacian/#","page":"4 p-Laplacian","title":"4 p-Laplacian","text":"This page was generated using Literate.jl.","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"EditURL = \"https://github.com/gridap/Tutorials/blob/master/src/inc_navier_stokes.jl\"","category":"page"},{"location":"pages/t008_inc_navier_stokes/#Tutorial-8:-Incompressible-Navier-Stokes-1","page":"8 Incompressible Navier-Stokes","title":"Tutorial 8: Incompressible Navier-Stokes","text":"","category":"section"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"(Image: ) (Image: )","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"In this tutorial, we will learn","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"How to solve nonlinear multi-field PDEs in Gridap\nHow to build FE spaces whose functions have zero mean value","category":"page"},{"location":"pages/t008_inc_navier_stokes/#Problem-statement-1","page":"8 Incompressible Navier-Stokes","title":"Problem statement","text":"","category":"section"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"The goal of this last tutorial is to solve a nonlinear multi-field PDE. As a model problem, we consider a well known benchmark in computational fluid dynamics, the lid-driven cavity for the incompressible Navier-Stokes equations. Formally, the PDE we want to solve is: find the velocity vector u and the pressure p such that","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"leftlbrace\nbeginaligned\n-Delta u + mathitRe (ucdot nabla) u + nabla p = 0 text in Omega\nnablacdot u = 0 text in Omega\nu = g text on partialOmega\nendaligned\nright","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"where the computational domain is the unit square Omega doteq (01)^d, d=2, mathitRe is the Reynolds number (here, we take mathitRe=10), and (w cdot nabla) u = (nabla u)^t w is the well known convection operator. In this example, the driving force is the Dirichlet boundary velocity g, which is a non-zero horizontal velocity with a value of g = (10)^t on the top side of the cavity, namely the boundary (01)times1, and g=0 elsewhere on partialOmega. Since we impose Dirichlet boundary conditions on the entire boundary partialOmega, the mean value of the pressure is constrained to zero in order have a well posed problem,","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"int_Omega q rm dOmega = 0","category":"page"},{"location":"pages/t008_inc_navier_stokes/#Numerical-Scheme-1","page":"8 Incompressible Navier-Stokes","title":"Numerical Scheme","text":"","category":"section"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"In order to approximate this problem we chose a formulation based on inf-sub stable Q_kP_k-1 elements with continuous velocities and discontinuous pressures (see, e.g., [1] for specific details). The interpolation spaces are defined as follows. The velocity interpolation space is","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"V doteq v in C^0(Omega)^d v_Tin Q_k(T)^d text for all TinmathcalT ","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"where T denotes an arbitrary cell of the FE mesh mathcalT, and Q_k(T) is the local polynomial space in cell T defined as the multi-variate polynomials in T of order less or equal to k in each spatial coordinate. Note that, this is the usual continuous vector-valued Lagrangian FE space of order k defined on a mesh of quadrilaterals or hexahedra. On the other hand, the space for the pressure is","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"beginaligned\nQ_0 doteq q in Q int_Omega q rm dOmega = 0 text with\nQ doteq q in L^2(Omega) q_Tin P_k-1(T) text for all TinmathcalT\nendaligned","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"where P_k-1(T) is the polynomial space of multi-variate polynomials in T of degree less or equal to k-1. Note that functions in Q_0 are strongly constrained to have zero mean value. This is achieved in the code by removing one degree of freedom from the (unconstrained) interpolation space Q and adding a constant to the computed pressure so that the resulting function has zero mean value.","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"The weak form associated to these interpolation spaces reads: find (up)in U_g times Q_0 such that r(up)(vq)=0 for all (vq)in V_0 times Q_0 where U_g and V_0 are the set of functions in V fulfilling the Dirichlet boundary condition g and 0 on partialOmega respectively. The weak residual r evaluated at a given pair (up) is the linear form defined as","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"r(up)(vq) doteq a((up)(vq))+ c(u)(v)","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"with","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"beginaligned\na((up)(vq)) doteq int_Omega nabla v cdot nabla u rm dOmega - int_Omega (nablacdot v) p rm dOmega + int_Omega q (nabla cdot u) rm dOmega\nc(u)(v) doteq int_Omega v \tcdot left( (ucdotnabla) u right) rm dOmega\nendaligned","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"Note that the bilinear form a is associated with the linear part of the PDE, whereas c is the contribution to the residual resulting from the convective term.","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"In order to solve this nonlinear weak equation with a Newton-Raphson method, one needs to compute the Jacobian associated with the residual r. In this case, the Jacobian j evaluated at a pair (up) is the bilinear form defined as","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"j(up)((delta u delta p)(vq)) doteq a((delta udelta p)(vq)) + rm dc(u)(delta uv)","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"where rm dc results from the linearization of the convective term, namely","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"rm dc(u)(delta uv) doteq int_Omega v cdot left( (ucdotnabla) delta u right) rm dOmega + int_Omega v cdot left( (delta ucdotnabla) u right) rm dOmega","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"The implementation of this numerical scheme is done in Gridap by combining the concepts previously seen for single-field nonlinear PDEs and linear multi-field problems.","category":"page"},{"location":"pages/t008_inc_navier_stokes/#Discrete-model-1","page":"8 Incompressible Navier-Stokes","title":"Discrete model","text":"","category":"section"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"We start with the discretization of the computational domain. We consider a 100times100 Cartesian mesh of the unit square.","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"using Gridap\nn = 100\ndomain = (0,1,0,1)\npartition = (n,n)\nmodel = CartesianDiscreteModel(domain,partition)","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"For convenience, we create two new boundary tags, namely \"diri1\" and \"diri0\", one for the top side of the square (where the velocity is non-zero), and another for the rest of the boundary (where the velocity is zero).","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"labels = get_face_labeling(model)\nadd_tag_from_tags!(labels,\"diri1\",[6,])\nadd_tag_from_tags!(labels,\"diri0\",[1,2,3,4,5,7,8])","category":"page"},{"location":"pages/t008_inc_navier_stokes/#FE-spaces-1","page":"8 Incompressible Navier-Stokes","title":"FE spaces","text":"","category":"section"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"For the velocities, we need to create a conventional vector-valued continuous Lagrangian FE space. In this example, we select a second order interpolation.","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"D = 2\norder = 2\nV = TestFESpace(\n reffe=:Lagrangian, conformity=:H1, valuetype=VectorValue{D,Float64},\n model=model, labels=labels, order=order, dirichlet_tags=[\"diri0\",\"diri1\"])","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"The interpolation space for the pressure is built as follows","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"Q = TestFESpace(\n reffe=:PLagrangian, conformity=:L2, valuetype=Float64,\n model=model, order=order-1, constraint=:zeromean)","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"With the options reffe=:PLagrangian, valuetype=Float64, and order=order-1, we select the local polynomial space P_k-1(T) on the cells TinmathcalT. With the symbol :PLagrangian we specifically chose a local Lagrangian interpolation of type \"P\". Using :Lagrangian, would lead to a local Lagrangian of type \"Q\" since this is the default for quadrilateral or hexahedral elements. On the other hand, constraint=:zeromean leads to a FE space, whose functions are constrained to have mean value equal to zero, which is just what we need for the pressure space. With these objects, we build the trial multi-field FE spaces","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"uD0 = VectorValue(0,0)\nuD1 = VectorValue(1,0)\nU = TrialFESpace(V,[uD0,uD1])\nP = TrialFESpace(Q)\n\nY = MultiFieldFESpace([V, Q])\nX = MultiFieldFESpace([U, P])","category":"page"},{"location":"pages/t008_inc_navier_stokes/#Nonlinear-weak-form-1","page":"8 Incompressible Navier-Stokes","title":"Nonlinear weak form","text":"","category":"section"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"The different terms of the nonlinear weak form for this example are defined following an approach similar to the one discussed for the p-Laplacian equation, but this time using the notation for multi-field problems.","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"const Re = 10.0\n@law conv(u,∇u) = Re*(∇u')⋅u\n@law dconv(du,∇du,u,∇u) = conv(u,∇du)+conv(du,∇u)\n\nfunction a(x,y)\n u, p = x\n v, q = y\n ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u)\nend\n\nc(u,v) = v⊙conv(u,∇(u))\ndc(u,du,v) = v⊙dconv(du,∇(du),u,∇(u))\n\nfunction res(x,y)\n u, p = x\n v, q = y\n a(x,y) + c(u,v)\nend\n\nfunction jac(x,dx,y)\n u, p = x\n v, q = y\n du, dp = dx\n a(dx,y)+ dc(u,du,v)\nend","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"With the functions res, and jac representing the weak residual and the Jacobian, we build the nonlinear FE problem:","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"trian = Triangulation(model)\ndegree = (order-1)*2\nquad = CellQuadrature(trian,degree)\nt_Ω = FETerm(res,jac,trian,quad)\nop = FEOperator(X,Y,t_Ω)","category":"page"},{"location":"pages/t008_inc_navier_stokes/#Nonlinear-solver-phase-1","page":"8 Incompressible Navier-Stokes","title":"Nonlinear solver phase","text":"","category":"section"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"To finally solve the problem, we consider the same nonlinear solver as previously considered for the p-Laplacian equation.","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"using LineSearches: BackTracking\nnls = NLSolver(\n show_trace=true, method=:newton, linesearch=BackTracking())\nsolver = FESolver(nls)","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"In this example, we solve the problem without providing an initial guess (a default one equal to zero will be generated internally)","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"uh, ph = solve(solver,op)","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"Finally, we write the results for visualization (see next figure).","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"writevtk(trian,\"ins-results\",cellfields=[\"uh\"=>uh,\"ph\"=>ph])","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"(Image: )","category":"page"},{"location":"pages/t008_inc_navier_stokes/#References-1","page":"8 Incompressible Navier-Stokes","title":"References","text":"","category":"section"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"[1] H. C. Elman, D. J. Silvester, and A. J. Wathen. Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics. Oxford University Press, 2005.","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"","category":"page"},{"location":"pages/t008_inc_navier_stokes/#","page":"8 Incompressible Navier-Stokes","title":"8 Incompressible Navier-Stokes","text":"This page was generated using Literate.jl.","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"EditURL = \"https://github.com/gridap/Tutorials/blob/master/src/poisson.jl\"","category":"page"},{"location":"pages/t001_poisson/#Tutorial-1:-Poisson-equation-1","page":"1 Poisson equation","title":"Tutorial 1: Poisson equation","text":"","category":"section"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"(Image: ) (Image: )","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"In this tutorial, we will learn","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"How to solve a simple PDE in Julia with Gridap\nHow to load a discrete model (aka a FE mesh) from a file\nHow to build a conforming Lagrangian FE space\nHow to define the different terms in a weak form\nHow to impose Dirichlet and Neumann boundary conditions\nHow to visualize results","category":"page"},{"location":"pages/t001_poisson/#Problem-statement-1","page":"1 Poisson equation","title":"Problem statement","text":"","category":"section"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"In this first tutorial, we provide an overview of a complete simulation pipeline in Gridap: from the construction of the FE mesh to the visualization of the computed results. To this end, we consider a simple model problem: the Poisson equation. We want to solve the Poisson equation on the 3D domain depicted in next figure with Dirichlet and Neumann boundary conditions. Dirichlet boundary conditions are applied on Gamma_rm D, being the outer sides of the prism (marked in red). Non-homogeneous Neumann conditions are applied to the internal boundaries Gamma_rm G, Gamma_rm Y, and Gamma_rm B (marked in green, yellow and blue respectively). And homogeneous Neumann boundary conditions are applied in Gamma_rm W, the remaining portion of the boundary (marked in white).","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"(Image: )","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"Formally, the problem to solve is: find the scalar field u such that","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"leftlbrace\nbeginaligned\n-Delta u = f textin Omega\nu = g texton Gamma_rm D\nnabla ucdot n = h texton Gamma_rm N\nendaligned\nright","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"being n the outwards unit normal vector to the Neumann boundary Gamma_rm N doteq Gamma_rm GcupGamma_rm YcupGamma_rm BcupGamma_rm W. In this example, we chose f(x) = 1, g(x) = 2, and h(x)=3 on Gamma_rm GcupGamma_rm YcupGamma_rm B and h(x)=0 on Gamma_rm W. The variable x is the position vector x=(x_1x_2x_3).","category":"page"},{"location":"pages/t001_poisson/#Numerical-scheme-1","page":"1 Poisson equation","title":"Numerical scheme","text":"","category":"section"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"To solve this PDE, we use a conventional Galerkin finite element (FE) method with conforming Lagrangian FE spaces (see, e.g., [1] for specific details on this formulation). The weak form associated with this formulation is: find uin U_g such that $ a(u,v) = b(v) $ for all vin V_0, where U_g and V_0 are the subset of functions in H^1(Omega) that fulfill the Dirichlet boundary condition g and 0 respectively. The bilinear and linear forms for this problems are","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":" a(uv) doteq int_Omega nabla v cdot nabla u rm dOmega quad b(v) doteq int_Omega v f rm dOmega + int_Gamma_rm N v h rm dGamma_rm N","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"The problem is solved numerically by approximating the spaces U_g and V_0 by their discrete counterparts associated with a FE mesh of the computational domain Omega. As we have anticipated, we consider standard conforming Lagrangian FE spaces for this purpose.","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"The implementation of this numerical scheme in Gridap is done in a user-friendly way thanks to the abstractions provided by the library. As it will be seen below, all the mathematical objects involved in the definition of the discrete weak problem have a correspondent representation in the code.","category":"page"},{"location":"pages/t001_poisson/#Setup-1","page":"1 Poisson equation","title":"Setup","text":"","category":"section"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"The step number 0 in order to solve the problem is to load the Gridap library in the code. If you have configured your Julia environment properly, it is simply done with the line:","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"using Gridap","category":"page"},{"location":"pages/t001_poisson/#Discrete-model-1","page":"1 Poisson equation","title":"Discrete model","text":"","category":"section"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"As in any FE simulation, we need a discretization of the computational domain (i.e., a FE mesh). All geometrical data needed for solving a FE problem is provided in Gridap by types inheriting from the abstract type DiscreteModel. In the following line, we build an instance of DiscreteModel by loading a json file.","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"model = DiscreteModelFromFile(\"../models/model.json\")","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"The file \"model.json\" is a regular json file that includes a set of fields that describe the discrete model. It was generated by using together the GMSH mesh generator and the GridapGmsh package. First, we generate a \"model.msh\" file with GMSH (which contains a FE mesh and information about user-defined physical boundaries in {GMSH} format). Then, this file is converted to the Gridap-compatible \"model.json\" file using the conversion tools available in the GridapGmsh package. See the documentation of the GridapGmsh for more information.","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"You can easily inspect the generated discrete model in Paraview by writing it in vtk format.","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"writevtk(model,\"model\")","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"The previous line generates four different files model_0.vtu, model_1.vtu, model_2.vtu, and model_3.vtu containing the vertices, edges, faces, and cells present in the discrete model. Moreover, you can easily inspect which boundaries are defined within the model.","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"For instance, if you want to see which faces of the model are on the boundary Gamma_rm B (i.e., the walls of the circular perforation), open the file model_2.vtu and chose coloring by the element field \"circle\". You should see that only the faces on the circular hole have a value different from zero (see next figure).","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"(Image: )","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"It is also possible to see which vertices are on the Dirichlet boundary Gamma_rm D. To do so, open the file model_0.vtu and chose coloring by the field \"sides\" (see next figure).","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"(Image: )","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"That is, the boundary Gamma_rm B (i.e., the walls of the circular hole) is called \"circle\" and the Dirichlet boundary Gamma_rm D is called \"sides\" in the model. In addition, the walls of the triangular hole Gamma_rm G and the walls of the square hole Gamma_rm Y are identified in the model with the names \"triangle\" and \"square\" respectively. You can easily check this by opening the corresponding file in Paraview.","category":"page"},{"location":"pages/t001_poisson/#FE-spaces-1","page":"1 Poisson equation","title":"FE spaces","text":"","category":"section"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"Once we have a discretization of the computational domain, the next step is to generate a discrete approximation of the finite element spaces V_0 and U_g (i.e. the test and trial FE spaces) of the problem. To do so, first, we are going to build a discretization of V_0 as the standard Conforming Lagrangian FE space (with zero boundary conditions) associated with the discretization of the computational domain. The approximation of the FE space V_0 is built as follows:","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"V0 = TestFESpace(\n reffe=:Lagrangian, order=1, valuetype=Float64,\n conformity=:H1, model=model, dirichlet_tags=\"sides\")","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"Here, we have used the TestFESpace constructor, which constructs a particular FE space (to be used as a test space) from a set of options described as key-word arguments. The with the options reffe=:Lagrangian, order=1, and valuetype=Float64, we define the local interpolation at the reference FE element. In this case, we select a scalar-valued, first order, Lagrangian interpolation. In particular, the value of the shape functions will be represented with 64-bit floating point numbers. With the key-word argument conformity we define the regularity of the interpolation at the boundaries of the cells in the mesh. Here, we use conformity=:H1, which means that the resulting interpolation space is a subset of H^1(Omega) (i.e., continuous shape functions). On the other hand, with the key-word argument model, we select the discrete model on top of which we want to construct the FE space. Finally, we pass the identifiers of the Dirichlet boundary via the dirichlet_tags argument. In this case, we mark as Dirichlet all objects of the discrete model identified with the \"sides\" tag. Since this is a test space, the corresponding shape functions vanishes at the Dirichlet boundary.","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"Once the space V_0 is discretized in the code, we proceed with the approximation of the trial space U_g.","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"g(x) = 2.0\nUg = TrialFESpace(V0,g)","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"To this end, we have used the TrialFESpace constructors. Note that we have passed a function representing the value of the Dirichlet boundary condition, when building the trial space.","category":"page"},{"location":"pages/t001_poisson/#Numerical-integration-1","page":"1 Poisson equation","title":"Numerical integration","text":"","category":"section"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"Once we have built the interpolation spaces, the next step is to set up the machinery to perform the integrals in the weak form numerically. Here, we need to compute integrals on the interior of the domain Omega and on the Neumann boundary Gamma_rm N. In both cases, we need two main ingredients. We need to define an integration mesh (i.e. a triangulation of the integration domain), plus a Gauss-like quadrature in each of the cells in the triangulation. In Gridap, integration meshes are represented by types inheriting from the abstract type Triangulation. For integrating on the domain Omega, we build the following triangulation and quadrature:","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"trian = Triangulation(model)\ndegree = 2\nquad = CellQuadrature(trian,degree)","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"Here, we build a triangulation from the cells of the model and define a quadrature of degree 2 in the cells of this triangulation. This is enough for integrating the corresponding terms of the weak form exactly for an interpolation of order 1.","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"On the other hand, we need a special type of triangulation, represented by the type\t BoundaryTriangulation, to integrate on the boundary. Essentially, a BoundaryTriangulation is a particular type of Triangulation that is aware of which cells in the model are touched by faces on the boundary. We build an instance of this type from the discrete model and the names used to identify the Neumann boundary as follows:","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"neumanntags = [\"circle\", \"triangle\", \"square\"]\nbtrian = BoundaryTriangulation(model,neumanntags)\nbquad = CellQuadrature(btrian,degree)","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"In addition, we have created a quadrature of degree 2 on top of the cells in the triangulation for the Neumann boundary.","category":"page"},{"location":"pages/t001_poisson/#Weak-form-1","page":"1 Poisson equation","title":"Weak form","text":"","category":"section"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"With all the ingredients presented so far, we are ready to define the weak form. This is done by means of types inheriting from the abstract type FETerm. In this tutorial, we will use the sub-types AffineFETerm and FESource. An AffineFETerm is a term that contributes both to the system matrix and the right-hand-side vector, whereas a FESource only contributes to the right hand side vector. Here, we use an AffineFETerm to represent all the terms in the weak form that are integrated over the interior of the domain Omega.","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"f(x) = 1.0\na(u,v) = ∇(v)⊙∇(u)\nb_Ω(v) = v*f\nt_Ω = AffineFETerm(a,b_Ω,trian,quad)","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"In the first argument of the AffineFETerm constructor, we pass a function that represents the integrand of the bilinear form a(cdotcdot). The second argument is a function that represents the integrand of the part of the linear form b(cdot) that is integrated over the domain Omega. The third argument is the Triangulation on which we want to perform the integration (in that case the integration mesh for Omega), and the last argument is the CellQuadrature needed to perform the integration numerically. Since the contribution of the Neumann boundary condition is integrated over a different domain, it cannot be included in the previous AffineFETerm. To account for it, we use a FESource:","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"h(x) = 3.0\nb_Γ(v) = v*h\nt_Γ = FESource(b_Γ,btrian,bquad)","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"In the first argument of the FESource constructor, we pass a function representing the integrand of the Neumann boundary condition. In the two last arguments we pass the triangulation and quadrature for the Neumann boundary.","category":"page"},{"location":"pages/t001_poisson/#FE-Problem-1","page":"1 Poisson equation","title":"FE Problem","text":"","category":"section"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"At this point, we can built the FE problem that, once solved, will provide the numerical solution we are looking for. A FE problem is represented in Gridap by types inheriting from the abstract type FEOperator (both for linear and nonlinear cases). Since we want to solve a linear problem, we use the concrete type AffineFEOperator, i.e., a problem represented by a matrix and a right hand side vector.","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"op = AffineFEOperator(Ug,V0,t_Ω,t_Γ)","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"Note that the AffineFEOperator object representing our FE problem is built from the test and trial FE spaces V0 and Ug, and the objects t_Ω and t_Γ representing the weak form.","category":"page"},{"location":"pages/t001_poisson/#Solver-phase-1","page":"1 Poisson equation","title":"Solver phase","text":"","category":"section"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"We have constructed a FE problem, the last step is to solve it. In Gridap, FE problems are solved with types inheriting from the abstract type FESolver. Since this is a linear problem, we use a LinearFESolver:","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"ls = LUSolver()\nsolver = LinearFESolver(ls)","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"LinearFESolver objects are built from a given algebraic linear solver. In this case, we use a LU factorization. Now we are ready to solve the FE problem with the FE solver as follows:","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"uh = solve(solver,op)","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"The solve function returns the computed numerical solution uh. This object is an instance of FEFunction, the type used to represent a function in a FE space. We can inspect the result by writing it into a vtk file:","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"writevtk(trian,\"results\",cellfields=[\"uh\"=>uh])","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"which will generate a file named results.vtu having a nodal field named \"uh\" containing the solution of our problem (see next figure).","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"(Image: )","category":"page"},{"location":"pages/t001_poisson/#References-1","page":"1 Poisson equation","title":"References","text":"","category":"section"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"[1] C. Johnson. Numerical Solution of Partial Differential Equations by the Finite Element Method. Dover Publications, 2009.","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"","category":"page"},{"location":"pages/t001_poisson/#","page":"1 Poisson equation","title":"1 Poisson equation","text":"This page was generated using Literate.jl.","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"EditURL = \"https://github.com/gridap/Tutorials/blob/master/src/elasticity.jl\"","category":"page"},{"location":"pages/t003_elasticity/#Tutorial-3:-Linear-elasticity-1","page":"3 Linear elasticity","title":"Tutorial 3: Linear elasticity","text":"","category":"section"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"(Image: ) (Image: )","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"In this tutorial, we will learn","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"How to approximate vector-valued problems\nHow to solve problems with complex constitutive laws\nHow to impose Dirichlet boundary conditions only in selected components\nHow to impose Dirichlet boundary conditions described by more than one function","category":"page"},{"location":"pages/t003_elasticity/#Problem-statement-1","page":"3 Linear elasticity","title":"Problem statement","text":"","category":"section"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"In this tutorial, we detail how to solve a linear elasticity problem defined on the 3D domain depicted in next figure.","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"(Image: )","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"We impose the following boundary conditions. All components of the displacement vector are constrained to zero on the surface Gamma_rm G, which is marked in green in the figure. On the other hand, the first component of the displacement vector is prescribed to the value deltadoteq 5mm on the surface Gamma_rm B, which is marked in blue. No body or surface forces are included in this example. Formally, the PDE to solve is","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"leftlbrace\nbeginaligned\n-cdotsigma(u) = 0 textin Omega\nu = 0 texton Gamma_rm G\nu_1 = delta texton Gamma_rm B\nsigma(u)cdot n = 0 texton Gamma_rm N\nendaligned\nright","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"The variable u stands for the unknown displacement vector, the vector n is the unit outward normal to the Neumann boundary Gamma_rm NdoteqpartialOmegasetminusleft(Gamma_rm BcupGamma_rm Gright) and sigma(u) is the stress tensor defined as","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"sigma(u) doteq lambda rm tr(varepsilon(u)) I +2 mu varepsilon(u)","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"where I is the 2nd order identity tensor, and lambda and mu are the Lamé parameters of the material. The operator varepsilon(u)doteqfrac12left(nabla u + (nabla u)^t right) is the symmetric gradient operator (i.e., the strain tensor). Here, we consider material parameters corresponding to aluminum with Young's modulus E=70cdot 10^9 Pa and Poisson's ratio nu=033. From these values, the Lamé parameters are obtained as lambda = (Enu)((1+nu)(1-2nu)) and mu=E(2(1+nu)).","category":"page"},{"location":"pages/t003_elasticity/#Numerical-scheme-1","page":"3 Linear elasticity","title":"Numerical scheme","text":"","category":"section"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"As in previous tutorial, we use a conventional Galerkin FE method with conforming Lagrangian FE spaces. For this formulation, the weak form is: find uin U such that $ a(u,v) = 0 $ for all vin V_0, where U is the subset of functions in VdoteqH^1(Omega)^3 that fulfill the Dirichlet boundary conditions of the problem, whereas V_0 are functions in V fulfilling v=0 on Gamma_rm G and v_1=0 on Gamma_rm B. The bilinear form of the problem is","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"a(uv)doteq int_Omega varepsilon(v) sigma(u) rm dOmega","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"The main differences with respect to previous tutorial is that we need to deal with a vector-valued problem, we need to impose different prescribed values on the Dirichlet boundary, and the integrand of the bilinear form a(cdotcdot) is more complex as it involves the symmetric gradient operator and the stress tensor. However, the implementation of this numerical scheme is still done in a user-friendly way since all these features can be easily accounted for with the abstractions in the library.","category":"page"},{"location":"pages/t003_elasticity/#Discrete-model-1","page":"3 Linear elasticity","title":"Discrete model","text":"","category":"section"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"We start by loading the discrete model from a file","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"using Gridap\nmodel = DiscreteModelFromFile(\"../models/solid.json\")","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"In order to inspect it, write the model to vtk","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"writevtk(model,\"model\")","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"and open the resulting files with Paraview. The boundaries Gamma_rm B and Gamma_rm G are identified with the names \"surface_1\" and \"surface_2\" respectively. For instance, if you visualize the faces of the model and color them by the field \"surface_2\" (see next figure), you will see that only the faces on Gamma_rm G have a value different from zero.","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"(Image: )","category":"page"},{"location":"pages/t003_elasticity/#Vector-valued-FE-space-1","page":"3 Linear elasticity","title":"Vector-valued FE space","text":"","category":"section"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"The next step is the construction of the FE space. Here, we need to build a vector-valued FE space, which is done as follows:","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"order = 1\n\nV0 = TestFESpace(\n reffe=:Lagrangian, order=order, valuetype=VectorValue{3,Float64},\n conformity=:H1, model=model, dirichlet_tags=[\"surface_1\",\"surface_2\"],\n dirichlet_masks=[(true,false,false), (true,true,true)])","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"As in previous tutorial, we construct a continuous Lagrangian interpolation of order 1. The vector-valued interpolation is selected via the option valuetype=VectorValue{3,Float64}, where we use the type VectorValue{3,Float64}, which is the way Gridap represents vectors of three Float64 components. We mark as Dirichlet the objects identified with the tags \"surface_1\" and \"surface_2\" using the dirichlet_tags argument. Finally, we chose which components of the displacement are actually constrained on the Dirichlet boundary via the dirichlet_masks argument. Note that we constrain only the first component on the boundary Gamma_rm B (identified as \"surface_1\"), whereas we constrain all components on Gamma_rm G (identified as \"surface_2\").","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"The construction of the trial space is slightly different in this case. The Dirichlet boundary conditions are described with two different functions, one for boundary Gamma_rm B and another one for Gamma_rm G. These functions can be defined as","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"g1(x) = VectorValue(0.005,0.0,0.0)\ng2(x) = VectorValue(0.0,0.0,0.0)","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"From functions g1 and g2, we define the trial space as follows:","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"U = TrialFESpace(V0,[g1,g2])","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"Note that the functions g1 and g2 are passed to the TrialFESpace constructor in the same order as the boundary identifiers are passed previously in the dirichlet_tags argument of the TestFESpace constructor.","category":"page"},{"location":"pages/t003_elasticity/#Constitutive-law-1","page":"3 Linear elasticity","title":"Constitutive law","text":"","category":"section"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"Once the FE spaces are defined, the next step is to define the weak form. In this example, the construction of the weak form requires more work than in previous tutorial since we need to account for the constitutive law that relates strain and stress. In this case, the integrand of the bilinear form of the problem is written in the code as follows:","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"a(u,v) = ε(v) ⊙ σ(ε(u))","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"The symmetric gradient operator is represented by the function ε provided by Gridap (also available as symmetric_gradient). However, function σ representing the stress tensor is not predefined in the library and it has to be defined ad-hoc by the user. The way function σ and other types of constitutive laws are defined in Gridap is by using the supplied macro @law:","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"using LinearAlgebra: tr\nconst E = 70.0e9\nconst ν = 0.33\nconst λ = (E*ν)/((1+ν)*(1-2*ν))\nconst μ = E/(2*(1+ν))\n@law σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"The macro @law is placed before a function definition. The arguments of the function annotated with the @law macro represent the values of different quantities at a generic integration point. In this example, the argument represents the strain tensor, from which the stress tensor is to be computed using the Lamé operator. Note that the implementation of function σ is very close to its mathematical definition. Under the hood, the @law macro adds an extra method to the annotated function. The newly generated method can be used as σ(ε(u)) in the definition of a bilinear form (as done above), or as σ(ε(uh)), in order to compute the stress tensor associated with a FEFunction object uh.","category":"page"},{"location":"pages/t003_elasticity/#Solution-of-the-FE-problem-1","page":"3 Linear elasticity","title":"Solution of the FE problem","text":"","category":"section"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"The remaining steps for solving the FE problem are essentially the same as in previous tutorial. We build the triangulation and quadrature for integrating in the volume, we define the terms in the weak form, and we define the FE problem. Finally, we solve it.","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"trian = Triangulation(model)\ndegree = 2*order\nquad = CellQuadrature(trian,degree)\nt_Ω = LinearFETerm(a,trian,quad)\nop = AffineFEOperator(U,V0,t_Ω)\nuh = solve(op)","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"Note that in the construction of the AffineFEOperator we have used a LinearFETerm instead of an AffineFETerm as it was done in previous tutorial. The LinearFETerm is a particular implementation of FETerm, which only leads to contributions to the system matrix (and not to the right hand side vector). This is what we need here since the body forces are zero. Note also that we do not have explicitly constructed a LinearFESolver. If a LinearFESolver is not passed to the solve function, a default solver is created and used internally.","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"Finally, we write the results to a file. Note that we also include the strain and stress tensors into the results file.","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"writevtk(trian,\"results\",cellfields=[\"uh\"=>uh,\"epsi\"=>ε(uh),\"sigma\"=>σ(ε(uh))])","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"It can be clearly observed (see next figure) that the surface Gamma_rm B is pulled in x_1-direction and that the solid deforms accordingly.","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"(Image: )","category":"page"},{"location":"pages/t003_elasticity/#Multi-material-problems-1","page":"3 Linear elasticity","title":"Multi-material problems","text":"","category":"section"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"We end this tutorial by extending previous code to deal with multi-material problems. Let us assume that the piece simulated before is now made of 2 different materials (see next figure). In particular, we assume that the volume depicted in dark green is made of aluminum, whereas the volume marked in purple is made of steel.","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"(Image: )","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"The two different material volumes are properly identified in the model we have previously loaded. To check this, inspect the model with Paraview (by writing it to vtk format as done before). Note that the volume made of aluminum is identified as \"material_1\", whereas the volume made of steel is identified as \"material_2\".","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"In order to build the constitutive law for the bi-material problem, we need a vector that contains information about the material each cell in the model is composed. This is achieved by these lines","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"using Gridap.Geometry\nlabels = get_face_labeling(model)\ndimension = 3\ntags = get_face_tag(labels,dimension)","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"Previous lines generate a vector, namely tags, whose length is the number of cells in the model and for each cell contains an integer that identifies the material of the cell. This is almost what we need. We also need to know which is the integer value associated with each material. E.g., the integer value associated with \"material_1\" (i.e. aluminum) is retrieved as","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"const alu_tag = get_tag_from_name(labels,\"material_1\")","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"Now, we know that cells whose corresponding value in the tags vector is alu_tag are made of aluminum, otherwise they are made of steel (since there are only two materials in this example).","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"At this point, we are ready to define the multi-material constitutive law. First, we define the material parameters for aluminum and steel respectively:","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"function lame_parameters(E,ν)\n λ = (E*ν)/((1+ν)*(1-2*ν))\n μ = E/(2*(1+ν))\n (λ, μ)\nend\n\nconst E_alu = 70.0e9\nconst ν_alu = 0.33\nconst (λ_alu,μ_alu) = lame_parameters(E_alu,ν_alu)\n\nconst E_steel = 200.0e9\nconst ν_steel = 0.33\nconst (λ_steel,μ_steel) = lame_parameters(E_steel,ν_steel)","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"Then, we define the function containing the constitutive law:","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"@law function σ_bimat(ε,tag)\n if tag == alu_tag\n return λ_alu*tr(ε)*one(ε) + 2*μ_alu*ε\n else\n return λ_steel*tr(ε)*one(ε) + 2*μ_steel*ε\n end\nend","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"Note that in this new version of the constitutive law, we have included a third argument that represents the integer value associated with a certain material. If the value corresponds to the one for aluminum (i.e., tag == alu_tag), then, we use the constitutive law for this material, otherwise, we use the law for steel.","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"Since we have constructed a new constitutive law, we need to re-define the bilinear form of the problem:","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"a(u,v) = ε(v) ⊙ σ_bimat(ε(u),tags)","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"In previous line, pay attention in the usage of the new constitutive law σ_bimat. Note that we have passed the vector tags containing the material identifiers in the last argument of the function`.","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"At this point, we can build the FE problem again and solve it","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"t_Ω = LinearFETerm(a,trian,quad)\nop = AffineFEOperator(U,V0,t_Ω)\nuh = solve(op)","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"Once the solution is computed, we can store the results in a file for visualization. Note that, we are including the stress tensor in the file (computed with the bi-material law).","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"writevtk(trian,\"results_bimat\",cellfields=\n [\"uh\"=>uh,\"epsi\"=>ε(uh),\"sigma\"=>σ_bimat(ε(uh),tags)])","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"Tutorial done!","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"","category":"page"},{"location":"pages/t003_elasticity/#","page":"3 Linear elasticity","title":"3 Linear elasticity","text":"This page was generated using Literate.jl.","category":"page"},{"location":"#Introduction-1","page":"Introduction","title":"Introduction","text":"","category":"section"},{"location":"#","page":"Introduction","title":"Introduction","text":"Welcome to the tutorial pages of the Gridap.jl project.","category":"page"},{"location":"#Contents-1","page":"Introduction","title":"Contents","text":"","category":"section"},{"location":"#","page":"Introduction","title":"Introduction","text":"Depth = 1","category":"page"},{"location":"#How-to-start-1","page":"Introduction","title":"How to start","text":"","category":"section"},{"location":"#","page":"Introduction","title":"Introduction","text":"There are different ways to use the tutorials:","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"[Recommended] Reading the html version of the tutorials. This is the recommended way if you want rapid access to the material with no setup steps. Simply click in one of the links in the Contents section.\n[Recommended] Running the Jupyter notebooks locally. A working installation of Julia in the system is required. See instructions in the How to run the notebooks locally section. This is the recommended way to follow the tutorials if you want to run the code and inspect the generated results with Paraview.\nRunning the notebook remotely via binder. In that case, go to the desired tutorial and click the icon (Image: ). No local installation of Julia needed.\nReading a non-interactive version of the notebook via nbviewer. In that case, go to the desired tutorial and click the icon (Image: )","category":"page"},{"location":"#How-to-run-the-notebooks-locally-1","page":"Introduction","title":"How to run the notebooks locally","text":"","category":"section"},{"location":"#","page":"Introduction","title":"Introduction","text":"Clone the repository","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"$ git clone https://github.com/gridap/Tutorials.git","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"Move into the folder and open a Julia REPL setting the current folder as the project environment. ","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"$ cd Tutorials\n$ julia --project=.\n _\n _ _ _(_)_ | Documentation: https://docs.julialang.org\n (_) | (_) (_) |\n _ _ _| |_ __ _ | Type \"?\" for help, \"]?\" for Pkg help.\n | | | | | | |/ _` | |\n | | |_| | | | (_| | | Version 1.1.0 (2019-01-21)\n _/ |\\__'_|_|_|\\__'_| | Official https://julialang.org/ release\n|__/ |\n\njulia> \n","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"Instantiate the environment. This will automatically download all required packages.","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"# Type ] to enter in pkg mode\n(Tutorials) pkg> instantiate","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"Build the notebooks","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"# Type Ctrl+C to get back to command mode\njulia> include(\"deps/build.jl\")","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"Open the notebooks","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"julia> using IJulia\njulia> notebook(dir=pwd())","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"This will open a browser window. Navigate to the notebooks folder and open the tutorial you want. Enjoy!","category":"page"},{"location":"#How-to-pull-the-latest-version-of-the-tutorials-1","page":"Introduction","title":"How to pull the latest version of the tutorials","text":"","category":"section"},{"location":"#","page":"Introduction","title":"Introduction","text":"If you have cloned the repository a while ago, you can update to the newest version with these steps.","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"Go to the Tutorials repo folder and git pull","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"$ git pull","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"Open Julia REPL","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"$ julia --project=.\n","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"and instantiate the environment and build the notebooks again","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"# Type ] to enter in pkg mode\n(Tutorials) pkg> instantiate\n\n# Type Ctrl+C to get back to command mode\njulia> include(\"deps/build.jl\")","category":"page"},{"location":"#","page":"Introduction","title":"Introduction","text":"Done!","category":"page"}]
}