Creating a New Element Type
The Element is the main part of a finite element program. It is here
that the program spend about 80% of it's time during calculation, but
when it comes to Explicit Finite Element codes, the task that the
element should perform is quite simple.
One can say the the task of a finite element is to do the following:
- Use the position of the nodes that the element is connected to,
to calculate the strain in the element.
- Use this strain and the material law that the element is supposed
to use, to figure out the resulting stress in the element.
- Use this stress and the geometry of the element to figure out
what the resulting reaction loads on the nodes will be, and apply them
to the nodes.
These three tasks are repeated over and over again in a simulation,
once for every timestep. If an element has more than one integration
point, the process has to be repeated once for every integration point
as well. This means that the calculation for an 8-integration point
solid element is 8 times longer per timestep, than for an 1-point solid.
This is also the reason why some fancy techniques such as hourglass
control have been developed to increase speed and allow use of fewer
integration points in elements.
If you want to add a new element, the challenge will be to find the
formulas you need to calculate the above tasks and then implement them.
The framework on how to do this is already within Impact and starts with
creating a new element type. It means creating a new sub-class to the
Element class. Examples of sub-classes are the Beam_2, Rod_2,
Shell_BT_4 and Solid_Iso_6 classes
Each subclass has the same set of methods (which are decided by the
Element class), but the coding is very different. We will investigate
each method later on in this chapter.
If you want to add a new element type, the easiest way is to pick a
type which can use the currently available material laws. At the time of
writing the material laws can handle:
- One dimensional stress
- Plane stress
- Three dimensional stress
The three dimensional stress is used in solid elements while the
plane stress is used in shell elements. One dimensional is used in the
Rod element.
You can read more about this in the material
chapter
The Element class is the mother class for all the
different element types. This class does not actually contain any code,
apart from a number of general methods which can be used by the
sub-classes. The rest of the methods are abstract methods. This means
that they act as a framework for the subclasses and if you don't include
each of the methods in your sub-class, you will get an error when
compiling.
The rest of Impact relies on that your new
sub-class has code for each of these methods. If you don't, you will not
be able to solve problems with your element.
Creating a sub-class is simple, but you should give it
a proper name. If you look at the names currently defined, they reflect:
When you create a new sub-class you must update the
list of elements in the getElementOfType_Fembic method and if
more file formats exists also in these related methods.
Each element sub-class also has a unique Type id,
stored in the Type field. The syntax for this parameter is a three digit
as follows:
First digit
The second digit is a group number (ex. rod has 0 and
Beam has 1)
The third digit is an index. If there are several
different rod types this is used to separate them
The Methods in the Element class are either
Abstract or General. General means that they are supposed to be used by
all the element sub-classes, but the coding resides in the element class
since it does not change from sub-class to sub-class.
The general methods are automatically inherited
and can be used by all sub-classes. Here follows a description:
This method is used by the reader class to
generate elements based on how they are described in the indata file.
The reason for putting this method here and not in the reader class is
once again to keep all the related work of adding an element type,
strictly within the element class. Therefore, the reader will use this
method to determine which element object to generate.
The _Fembic part of the name is related to the
fact that the Fembic type indata format is assumed. If a new indata file
format is added, i.e. a new reader sub-class is written, an extra
method will be required here, which is tailored to that file format.
Some elements use local coordinate systems to make
calculations simpler. Examples include shell elements. This method can
be used to calculate the transformation matrix that transform from local
to global coordinates system. The input are three nodal coordinates
which form the coordinate system axes. See the code for more detailed
explanation
This method searches the database and returns a
handle to the material object with a specific name. It is used by the
parse method in the element sub-classes to interpret element indata.
This method searches the database and returns a
handle to the node object with a specific number. It is used by the
parse method in the element sub-classes to interpret element indata.
This is quite a useful method. When you have a
string with a range of numbers that are separated by commas, this method
picks the n:th number in this string, converts it into an integer and
returns it. It can be useful when reading indata and is used by the
element sub-classes in their parse method.
Sometimes it can be useful to know if an element
has been asked or not. The parameter 'processed' can be set and read to
determine just that. This method returns the value of that parameter and
since it is inherited to all sub-elements, it is placed here in the
mother element class.
Allows you to set or clear that parameter
mentioned above.
The Abstract methods are empty in coding in the
element class. This is because every sub-class has a different
implementation of these methods. As a programmer that wishes to add a
new element sub-class, these methods are the ones that he/she has to
write. Here follows a description of each method
This method calculates the elements contribution
regarding mass and rotational inertia to the nodes and is called at the
initialisation of the problem. Since this program uses lumped mass, all
the mass are concentrated to the nodes, and so are the inertia of
rotation as well. The mass is usually quite simple to calculate. In the
example of a rod, it is half of the total mass of the rod in each node.
For a shell element it could be slightly more
trickier. The difficult part comes when the rotational inertia is to be
calculated since this is needs to be calculated in three dimensions and
then transformed to the global xyz coordinate system before adding it to
the node.
Roughly, the procedure is as follows:
- Calculate element mass distribution and add to
node
- Calculate element inertia
- Transform to global directions
- Add to node inertia matrix
The contact forces are calculated by this method.
Contact forces are forces that are generated onto nodes because they are
about to contact an element. To prevent the node to just pass through
the element, an artificial force is calculated and put on to the node
(with the opposite force put onto the element).
This is done in the following way:
- Let the element check it's segment(s) or faces
against any neighbouring nodes.
- calculate distance from segment to node.
- If distance is smaller than tolerance then calculate a reaction
force as a function of the distance.
- If a reaction force has been calculated, then add this force to
the node and the opposite force to the elements nodes.
- Continue this loop until all the nodes within the contact
tolerance from the element has been checked.
It is worth noting that there are loads of
different contact algorithms out there. Each of them with different
benefits. The good thing with this program is that each node knows about
it's closest neighbour node! This can be used to limit the amount of
nodes you need to check for your contact algorithm.
After each timestep, the nodes will update
themselves to see if their closes neighbour has changed.
This method computes the contribution from gravity
and external forces on the elements. The gravity and external forces on
the elements are transformed into nodal loads which are then added to
the nodes connected to the element.
This method calculates and adds the element
internal forces to the nodes. Please note that by definition, the
internal forces are subtracted from the external loads. Thus, the
"addition" is negative.
The coding can be quite extensive in this method
for some elements. In principle, the approach is however more or less
the same. If we look at the rod, the coding is quite simple:
Start by calculating the local force using stress
in rod and cross sectional area
force.set(0, 0,stress.get(0,0)*cross_section_area);
force.set(1, 0, 0);
force.set(2, 0, 0);
Now, continue by transforming this force to global
coordinates using a transformation matrix calculated in another element
method.
global_force = local_coordinate_system.times(force).copy();
Finally, add this force contribution to the nodes.
Remember to use negative addition since this is an internal force (a
reaction force).
node1.addForce(global_force.times(1));
node2.addForce(global_force.times(-1));
The strain of the element is calculated here. It
is a matter of transforming the positions of the nodes that the element
is connected to, into a strain matrix for the element.
In addition, the strain increment is also
calculated (that is the difference in strain between now and the strain
in the previous timestep)
Calculating the strain for the rod element is
quite simple. Start by calculating the length of the rod.
xpos1 = node1.getX_pos();
ypos1 = node1.getY_pos();
zpos1 = node1.getZ_pos();
//
xpos2 = node2.getX_pos();
ypos2 = node2.getY_pos();
zpos2 = node2.getZ_pos();
new_length = java.lang.Math.sqrt((xpos2 - xpos1) * (xpos2 - xpos1) +
(ypos2 - ypos1) * (ypos2 - ypos1) + (zpos2 - zpos1) * (zpos2 - zpos1));
The strain increment can now be calculated. Note
that his is true strain.
dstrain.set(0,0,Math.log(1+(new_length-initial_length)/initial_length) -
strain.get(0,0));
Finally, calculate the cross section area of the
rod, assuming incompressible material.
cross_section_area =
initial_cross_section_area*initial_length/new_length;
The strain matrix will be updated when the
material law is used to calculate the stress.
This method calculates the stress matrix of the
element out of a strain matrix and a strain increment matrix. The main
work here is actually done inside the material law of which the element
uses. Which law to use is assigned to the element in the initialisation
stage.
This is the code in the rod element. The stress
matrix is supplied as an indata, but is actually updated inside the
material law. Note that for the rod element, the one dimensional law is
used. This has to be implemented for each new material law that is added.
material.calculateStressOneDimensional(strain,dstrain,stress);
For each element there is a lower limit on the
critical timestep. This timestep is decided by the shortest travelling
distance a wave can take through the element on one step in time.
For a rod element, this is the length of the
element, divided by the wave-speed of the rod material. For a hexahedron
(solid) element, the critical length is the shortest side of the element
and thus, the timestep is this critical length divided by the wave-speed
of the material that the solid is made out of.
The material law supplies the wave-speed. This
method has input and if the element has a smaller timestep than the
input value, the element should return it's own timestep value.
This method is called by the main routine in
smack. The program is asking all elements for their timestep before it
takes it's next step. The size of this step is set by the element that
requires the smallest timestep. It sets the stability of the solution.
This method simply returns the number of the
element. This number is the one that was defined for the element in the
indata file.
A certain element type can have several different
number of integration points. An element with many interaction points has
greater accuracy and stability but takes longer time to calculate.
As an example, the solid element can either have
one interaction point or eight. The one integration point version is not
stable since so called hourglassing can occur. There are stabilisation
routines available that can fix this but at the time of writing, there
is none implemented.
Each element is assigned a type number (see
getElementOfType_Fembic above). This method returns this number.
As described before, some of the interpretation of
indata files has been put into the element itself. This method is called
when a Fembic indata file is read and the indata for the element needs
to be read and interpreted.
The beauty with this approach is that the element
can read the indata string and set all its variables by itself. The
variables can then be kept private and encapsulated inside the element.
Another benefit is that the programmer only needs
to work in this class when he/she adds a new element. There is no need
to fiddle around els where in impact.
When a new file format is added, a new additional
method is required for this file format and the programmer needs to add
this new method to all element types in impact.
This method is the opposite of the parse_Fembic
method described above. It reads all the parameters in the element and
creates an output string which is then printed into the outdata file.
This very method prints the outdata from the element in GID format. GID
is a pre- and postprocessor which is free and a good complement to
Impact.
This method is called every time the writer object
wants to print something, but there could be different things from time
to time like header data or bulk data. A control parameter is supplied
with the call to signal what is to be printed.
This is the big initialisation method for the
element and it is called once at the beginning of the solution. Examples
of what this method should do is:
If we look at how the material object is
initialised, here is example code:
try {
material = (Material) material.clone();
} catch (CloneNotSupportedException e) {
System.err.println("Object cannot clone");
}
The material object is told to clone itself and
the handle is reinitialised
This method simply assigns a number to the element.
Some elements has a local coordinate system. The
system is really a transformation matrix between the local and global
coordinates system and needs to be updated between each timestep since
the element deforms and the coordinate system deforms with it.
Here is a code example. The method
'calculateLocalBaseVectors' in the Element mother-class is used and is
supplied the element node coordinates.
local_coordinate_system = new
Jama.Matrix(super.calculateLocalBaseVectors(node1.getX_pos(),node1.getY_pos(),no
de1.getZ_pos(),node2.getX_pos(),node2.getY_pos(),node2.getZ_pos(),node2.getX_pos
()+1,node2.getY_pos()+1,node2.getZ_pos()+1).getArray());
Have a good look at the coding of the Rod_2 and the
Solid_Iso_6 element. The code is full of comments. Good luck!