Merge pull request #2 from JuliaFEM/release/v0.1.0

Implemented distributed coupling element

examples/example1.jl

@@ -0,0 +1,122 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMCoupling.jl/blob/master/LICENSE
+
+using JuliaFEM
+using FEMBase
+using Logging
+Logging.configure(level=DEBUG)
+
+type Coupling <: BoundaryProblem
+    reference_node :: Element{Poi1}
+end
+
+function Coupling()
+    reference_node = Element(Poi1, Int64[])
+    return Coupling(reference_node)
+end
+
+d = 3.0
+X = Dict(1 => [0.0, 0.0],
+         2 => [d, 0.0],
+         3 => [d, d],
+         4 => [0.0, d],
+         5 => [3*d, 1/2*d])
+element1 = Element(Quad4, [1, 2, 3, 4])
+update!(element1, "geometry", X)
+update!(element1, "youngs modulus", 288.0)
+update!(element1, "poissons ratio", 1/3)
+
+problem = Problem(Elasticity, "test problem", 2)
+problem.properties.formulation = :plane_strain
+add_elements!(problem, [element1])
+
+bc = Problem(Dirichlet, "fixed", 2, "displacement")
+element2 = Element(Seg2, [4, 1])
+update!(element2, "geometry", X)
+update!(element2, "displacement 1", 0.0)
+update!(element2, "displacement 2", 0.0)
+add_elements!(bc, [element2])
+
+"""
+    add_coupling_nodes!(problem, elements)
+
+Add new coupling nodes into the problem. Nodes must be defined with Poi1
+elements.
+"""
+function add_coupling_nodes!(problem, elements)
+    for element in elements
+        push!(problem.elements, element)
+    end
+    return nothing
+end
+
+"""
+    add_reference_node!(problem, element)
+
+Add new reference node into the problem. Node must be defined with a Poi1
+element. If reference node is already defined, it will be replaced with new one.
+"""
+function add_reference_node!(problem, element)
+    problem.properties.reference_node = element
+    return nothing
+end
+
+function FEMBase.assemble_elements!(problem::Problem{Coupling},
+                                    assembly::Assembly,
+                                    elements::Vector{Element{Poi1}}, time)
+
+    ref_node = problem.properties.reference_node
+    if length(get_connectivity(ref_node)) == 0
+        error("Reference node node defined. Define reference node using add_reference_node!")
+    end
+    ref_node_id = first(get_connectivity(ref_node))
+    info("Reference node id = $ref_node_id")
+    X_r = first(ref_node("geometry", time))
+    info("Reference node geometry = $X_r")
+    fe = zeros(2)
+    for coupling_node in elements
+        fill!(fe, 0.0)
+        coupling_node_id = first(get_connectivity(coupling_node))
+        X_n = first(coupling_node("geometry", time))
+        info("Coupling node id = $coupling_node_id, geometry = $X_n")
+        gdofs = get_gdofs(problem, coupling_node)
+        info("gdofs = $gdofs")
+        if haskey(ref_node, "point moment 3")
+            M = ref_node("point moment 3")
+            info("calculate point moment")
+            # fe += ...
+        end
+        if haskey(ref_node, "point load 1")
+            M = ref_node("point load 1")
+            info("calculate point moment")
+            # fe += ...
+        end
+        if haskey(ref_node, "point load 2")
+            M = ref_node("point load 2")
+            info("calculate point moment")
+            # fe += ...
+        end
+        add!(assembly.f, gdofs, fe)
+    end
+end
+
+coupling = Problem(Coupling, "test", 2, "displacement")
+# "slave nodes"
+element3 = Element(Poi1, [2])
+element4 = Element(Poi1, [3])
+update!([element3, element4], "geometry", X)
+element5 = Element(Poi1, [5])
+update!(element5, "geometry", X)
+update!(element5, "point moment 3", 10.0/d)
+add_coupling_nodes!(coupling, [element3, element4])
+add_reference_node!(coupling, element5)
+
+step = Analysis(Nonlinear)
+#xdmf = Xdmf("example1_results"; overwrite=true)
+#add_results_writer!(step, xdmf)
+add_problems!(step, [problem, bc, coupling])
+run!(step)
+#close(xdmf.hdf)
+time = 0.0
+u = element1("displacement", time)
+println("displacement: $u")

examples/example_dcoupling_cylindrical_beam.jl

@@ -0,0 +1,129 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMCoupling.jl/blob/master/LICENSE
+
+# # Cylinder in torsion using distributed coupling
+
+# Distributed coupling can be used to distribute point moments as point forces
+# to nodes which dont have rotational degrees of freedom.
+
+# The model is a 3d cylinder, shown in picture.
+
+# ![](example_dcoupling_cylindrical_beam/example_dcoupling_cylindrical_beam.png)
+
+# Distributed coupling is used to distribute a torque load to the end of the
+# cylinder. The other end is fixed. Node sets for coupling nodes, reference
+# node and fixed face were made in the ABAQUS input file. Reference node is
+# located in the centrum of the coupling nodes.
+
+using JuliaFEM
+using JuliaFEM: add_elements!, Problem
+using JuliaFEM.Preprocess
+using JuliaFEM.Abaqus: create_surface_elements
+using FEMBase
+using FEMCoupling
+using FEMCoupling: add_reference_node!, add_coupling_nodes!
+
+# reading mesh from ABAQUS input file
+
+datadir = Pkg.dir("FEMCoupling", "examples", "example_dcoupling_cylindrical_beam")
+mesh = abaqus_read_mesh(joinpath(datadir, "example_dcoupling_cylindrical_beam.inp"))
+println("Number of nodes in a model: ", length(mesh.nodes))
+
+# # Elements
+# Creating elements for the body. Mesh and element types are defined in
+# the ABAQUS input file. The cylinder body is named "Body1" in the input file.
+
+cylinder_body = create_elements(mesh,"Body1")
+
+# Updating values for the elements.
+
+update!(cylinder_body, "youngs modulus", 210e3)
+update!(cylinder_body, "poissons ratio", 0.3)
+
+# Creating an elasticity problem and adding the elements to it.
+
+cylinder_problem = Problem(Elasticity,"cylinder_problem",3)
+add_elements!(cylinder_problem, cylinder_body)
+
+# # Boundary conditions
+# Creating Poi1-type elements as boundary condition elements to nodes of the
+# node set Fixed_face_set.
+
+bc_elements = [Element(Poi1, [j]) for j in mesh.node_sets[:Fixed_face_set]]
+
+# Updating geometry for the bc elements
+
+update!(bc_elements, "geometry", mesh.nodes)
+
+# Fixing all displacements for the bc elements.
+
+for i=1:3
+    update!(bc_elements, "displacement $i", 0.0)
+end
+
+# Creating a bc problem and adding the bc elements to it.
+
+bc = Problem(Dirichlet, "fixed", 3, "displacement")
+add_elements!(bc, bc_elements)
+
+# # Distributed coupling
+# Creating Poi1 elements to nodes in coupling nodes set.
+
+coupling_nodes = [Element(Poi1, [j]) for j in mesh.node_sets[:Coupling_nodes_set]]
+
+# Updating geometry for the coupling nodes.
+
+update!(coupling_nodes, "geometry", mesh.nodes)
+
+# Creating Poi1 element for the reference node.
+
+reference_node_id = collect(mesh.node_sets[:ref_node_set])
+reference_node = Element(Poi1, reference_node_id)
+
+# Updating geometry and applying a point moment for the reference node.
+
+update!(reference_node, "geometry", mesh.nodes)
+update!(reference_node, "point moment 3", 1500.0)
+
+# Creating a coupling problem and adding coupling nodes and reference nodes to
+# it.
+
+coupling = Problem(Coupling, "cylind", 3, "displacement")
+add_coupling_nodes!(coupling, coupling_nodes)
+add_reference_node!(coupling, reference_node)
+
+# # Analysis
+# Creating a step and running the analysis. The cylinder_problem contains
+# information about the body, its elements and their values. The bc contains
+# information about boundary conditions and coupling contains information
+# about distributed coupling.
+
+step = Analysis(Nonlinear)
+add_problems!(step, [cylinder_problem, bc, coupling])
+run!(step)
+
+# # Results
+
+# Comparing calculated results with ABAQUS results. The node set
+# circlenodes_set contains nodes which are on the outer face radius.
+# These circle nodes should have the maximum displacement magnitude (norm(u)).
+
+node_on_circle = first(mesh.node_sets[:circlenodes_set])
+
+# declaring displacements at time 0.0 to variable u
+
+time=0.0
+u = cylinder_problem("displacement", time)[node_on_circle]
+u_mag = norm(u)
+
+# Making a testset.
+
+using FEMBase.Test
+@testset "displacement magnitude" begin
+u_mag_expected=6.306e-4
+@test isapprox(u_mag, u_mag_expected, rtol=1e-3)
+end
+
+# Printing node ids
+println("reference node id = $(reference_node_id[1])")
+println("node on circle id = $node_on_circle")