confusion about vectors sizes

[mahdimasmoudi]

I am trying to undestand few steps from the code and I need help understanding how results are derived.

I start with 45 nodes of a simple cube. I get 100 elements from this mesh.

A=get_matrix gives a 1785 x 1785 matrix b=get_vector gives a 1785 elements vector. I suppose I am going to solve the system AU=b where U is the displacement I am solving for.

The confusion comes out of the fact that U should have the dimension: number of nodes * number of DoFs which is 135 in this case and not 1785.

When I reach the last step to store the results in a vtk file it turns out that the vector stored is a 400 elements vector which is even more confusing.

I appreciate any help in this regard.

[JordiManyer]

@mahdimasmoudi I would probably need to see the code.

[mahdimasmoudi]

I am following a similar code of this tutorial. I am just changing the input msh file to a simple cube. https://gridap.github.io/Tutorials/dev/pages/t003_elasticity/#Tutorial-3:-Linear-elasticity-1

using GridapGmsh
using Gridap

model = GmshDiscreteModel("C:/Users/masmoud2/Desktop/parameters identification/cube small.msh")

writevtk(model,"model")

order = 3

reffe = ReferenceFE(lagrangian,VectorValue{3,Float64},order)
V0 = TestFESpace(model,reffe;
  conformity=:H1,
  dirichlet_tags=["fixed end "],
  dirichlet_masks=[(true,true,true)]);

g(x) = VectorValue(0.0,0.0,0.0)
U = TrialFESpace(V0,[g])

E = 100.0e9
ν = 0.4
λ = (E*ν)/((1+ν)*(1-2*ν))
μ = E/(2*(1+ν))
σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε

degree = 2*order
Ω = Triangulation(model)

dΩ = Measure(Ω,degree)

forectags = ["forced side"]
Γ = BoundaryTriangulation(model,tags=forectags)
dΓ = Measure(Γ,degree)

f(x) = VectorValue( 0.0, 1.0e9,0.0) # external force
a(u,w) = ∫( ε(w) ⊙ (σ∘ε(u)) )*dΩ
l(v) =  ∫(v⋅f)*dΓ
op = AffineFEOperator(a,l,U,V0)
op.op.vector  # this has a size of 1785 
op.op.matrix

uh = Gridap.solve(op)

writevtk(Ω,"results5",cellfields=["uh"=>uh])

[JordiManyer]

Is the mesh composed of cubes or tetrahedra?

[JordiManyer]

Let's assume it's cubes, i.e you have N_cells = 100 cuboid cells.

In that case, each reference element has N_nodes_x_element = (1+order)^D nodes with D components each. In your particular case, 64 nodes and 192 dofs per cell. Some of them are indeed shared between cells (the ones in the nodes, edges and faces) but some of them are not (the interior ones).

On top of this, you have dirichlet boundary conditions. That means that all dofs on the boundary will be eliminated and their contributions with go to your rhs. So those get subtracted from your total matrix size.

When writing to vtk, you are not using the option to store higher-dimensional data. So your solution is interpolated back into first order. Probably the dirichlet dofs are added back, since VTK needs complete nodal data.

So it's quite a bit more complicated than you think. If you want the simplest case, I would recommend a cartesian mesh with first order 2D scalar elements without any boundary conditions. You'll probably get easier counts.