Simple linear example

[sebastianpech]

I wanted to run a simple linear elastic beam example, however, I didn't manage to get the analysis to run. I used test/test_beam3d_3_beams.jl, test/test_supports.jl and the modal analysis of the frame in JuliaFEM as templates. This is my code:

using FEMBeam, JuliaFEM

# Define nodal coordinates
L = 10
N = 10
X = Dict(n => [0.0 + (n-1)*L/(N-1), 0.0, 0.0] for n in 1:N)

# Create beams
beams = Element.(Seg2, [[n, n+1] for n in 1:N-1])

# Define beam properties
update!(beams, "geometry", X)
update!(beams, "youngs modulus", 100.0e6)
update!(beams, "shear modulus", 80.0e6)
update!(beams, "cross-section area", 2.5e-2)
update!(beams, "torsional moment of inertia 1", 1.0e-4)
update!(beams, "torsional moment of inertia 2", 1.5e-4)
update!(beams, "polar moment of inertia", 3.0e-5)
update!(beams, "normal", [0.0, 0.0, -1.0])

# Define BC. Clamp node 1
fix = Element(Poi1, [1])
for i=1:3
    global fix
    update!(fix, "fixed displacement $i", 0.0)
    update!(fix, "fixed rotation $i", 0.0)
end

# Define BC. Add vertical point load at Node N
load = Element(Poi1, [N])
update!(load, "point force 3",  1.0e3)

# Define Problem
problem = Problem(Beam, "example 1", 6)
add_elements!(problem, beams)
add_elements!(problem, fix)
add_elements!(problem, load)

# Define Analysis
analysis = Analysis(Linear, problem)

# Define output and run analysis
xdmf = Xdmf("model_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
push!(problem.postprocess_fields, "stress")
run!(analysis)
close(xdmf)

The script failes with

 Warning: No rows in C2, forget to set Dirichlet boundary conditions to model?
└ @ JuliaFEM ~/.julia/packages/JuliaFEM/DbhTT/src/solvers.jl:197
ERROR: LoadError: ArgumentError: the number of rows or columns of the matrix must be greater than zero

Can you please check what I'm missing.

[sebastianpech]

I also tried the version with the extra problem definition for the dirichlet BC, but I get the same error message.

[ahojukka5]

For me everything looks good. It looks that Dirichlet boundary conditions are missing. Here is the code for handling boundary conditions. What about if you do

assemble!(problem, 0.0)
show(sparse(problem.assembly.C2))

There should be ones in diagonal for fixed dofs and it should be same than problem.assembly.C1.

[sebastianpech]

There should be ones in diagonal for fixed dofs and it should be same than problem.assembly.C1.

That works as expected:

julia> show(sparse(problem.assembly.C2))

  [1, 1]  =  1.0
  [2, 2]  =  1.0
  [3, 3]  =  1.0
  [4, 4]  =  1.0
  [5, 5]  =  1.0
  [6, 6]  =  1.0
julia> show(sparse(problem.assembly.C1))

  [1, 1]  =  1.0
  [2, 2]  =  1.0
  [3, 3]  =  1.0
  [4, 4]  =  1.0
  [5, 5]  =  1.0
  [6, 6]  =  1.0

[sebastianpech]

I debugged the get_boundary_assembly function and for the code above it returns C1 and C2 beeing 0. I think the loop in get_boundary_assembly is not running because JuliaFEM.get_boundary_problems(analysis) returns an empty array due to this.

[ahojukka5]

Hmm yes. I'm not sure why we even have FieldProblem and BoundaryProblem types separately. This is not the best kind of design. In earlier stage, we thought that we could have separate boundary conditions assembly which can be used in all problem types, but it's not too much of boilerplate to have boundary assembly done in each problem type separately.

What about, if you do in the script, before assembling,

import FEMBase
FEMBase.is_field_problem(::Problem{Beam}) = true
FEMBase.is_boundary_problem(::Problem{Beam}) = true
@assert is_field_problem(problem)
@assert is_boundary_problem(problem)

does this solve the issue?

[sebastianpech]

Partly. I also had to define FEMBase.get_parent_field_name(p::Problem{Beam}) = p.parent_field_name. However, the post processing methods are not defined

ERROR: MethodError: no method matching postprocess!(::Problem{Beam}, ::Float64, ::Type{Val{:stress}})
Closest candidates are:
  postprocess!(::Problem{Elasticity}, ::Float64, ::Type{Val{:stress}}) at /Users/sebastianpech/.julia/packages/JuliaFEM/DbhTT/src/problems_elasticity.jl:468
  postprocess!(::Problem{Elasticity}, ::Float64, ::Type{Val{:strain}}) at /Users/sebastianpech/.julia/packages/JuliaFEM/DbhTT/src/problems_elasticity.jl:462
  postprocess!(::Problem{Dirichlet}, ::Float64, ::Type{Val{Symbol("reaction force")}}) at /Users/sebastianpech/.julia/packages/JuliaFEM/DbhTT/src/problems_dirichlet.jl:157

Are there any defined for Beams?

[ahojukka5]

No, there is not. I don't see any postprocessing code here. Just a basic implementation of Euler-Bernoulli beams. I wrote some implementation details here. It would not be hard to implement calculating stress, though.

[sebastianpech]

Sure. I removed them for now and I get plausible results. Thanks for the help

[ahojukka5]

So basically we should add

FEMBase.is_field_problem(::Problem{Beam}) = true
FEMBase.is_boundary_problem(::Problem{Beam}) = true

to somewhere to FEMBase.jl code to solve this issue.

[sebastianpech]

You also need

FEMBase.get_parent_field_name(p::Problem{Beam}) = p.parent_field_name

because it's just defined for a BoundaryProblem here.

[sebastianpech]

Should get_parent_field_name be implemented for non-boundary problems as well?

[ahojukka5]

I think no. The idea behind this "parent field" was that when the primary field of Dirichlet boundary condition is "lambda" as we are basically building a mixed problem, we needed somehow to transfer the information what is the name of the actual field considered to be solved. So we can use the same Dirichlet for both Elasticity and Heat. But I think this is an over-complicated solution and benefits are minimal. We should have only one kind of problem instead of FieldProblem and BoundaryProblem. All problems should also take care of building the boundary conditions like is done in this package, it seems to be a better way.

[ahojukka5]

We use get_parent_field_name only in a couple of places.

• https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/src/problems_dirichlet.jl#L67
• https://github.com/JuliaFEM/FEMBase.jl/blob/master/src/problems.jl#L352

And in some contact algorithms. In those, it's important to know is the contact formulation done for displacement or temperature, but I guess there is better ways to organize that.

[sebastianpech]

What do you mean by

All problems should also take of building the boundary conditions like is done in this package, it seems to be a better way.

With

update!(fix, "fixed displacement 1", 0.0)
update!(fix, "fixed rotation 1", 0.0)

instead of an extra Dirichlet problem?

[ahojukka5]

I think each package should deal with its own boundary conditions like is done in this package. That makes it much easier to handle all kinds of special cases. For example, for beams, we could have slanted supports and things like that. And the right place for implementation of those is FEMBeams.jl.

[sebastianpech]

Ah, ok

[sebastianpech]

I just noticed that all dofs in the xmf file are collected in a single point array (so displacement also contains the rotations). This causes a problem when viewing the results because for warping paraview needs a three point array. Is it possible to separate the dofs?

[ahojukka5]

I did a dirty workaround when storing eigenmodes. Maybe something similar dirty workaround to here would do the trick. Paraview cannot use or visualize rotation dofs. Another option would be to study Paraview documentation, probably it is possible to apply some kind of filter to the dataset before visualization.

[sebastianpech]

Ok, thanks for your help

[sebastianpech]

I ended up writing my own script using Plots.jl.