How are Dirichlet boundary conditions implemented in FEAP

[festein]

Hello,
I currently use FEAP to check self-developed FEM implementations whether they are correct or not. While doing so, I wondered how exactly FEAP implements Dirichlet boundary conditions. It is important for me since I want to compare some internal element values between implementations, so I need to make sure, that my own implementations handle the Dirichlet boundary conditions the same way FEAP does.

Generally, I am familiar with 3 types of implementations
1. Adding the Dirichlet boundary conditions to the displacement vector when developing the element stiffness matrices and residuals. Then the respective degrees of freedom are removed from the global system of equations before the solve step.
2. Disregarding the Dirichlet boundary conditions in the assembly and then adding them to the RHS of the global system by multiplying their values with the respective columns of the global stiffness matrix. Furthermore, the respective degrees of freedom are removed and the system is solved.
3. Lagrange parameters.
I experimented a bit with FEAP and if I use PRINT,FORM I can see that the size of my residual is as big as the reduced system (global degrees of freedom - Dirichlet boundary degrees of freedom). Therefore option 1 or 2 needs to be used. But how is it done exactly and when (in which batch file step) does it happen? Also, is there a specific reason why the certain procedure is chosen?
Kind regards,
Felix

[JStorm]

FEAP makes use of option number 2.
The element stiffness is directly assembled into the reduced system of equations and the RHS vector is modified by the disregarded stiffness term and the boundary condition value.

[Prof. R.L. Taylor]

FEAP uses Option 1:  After leaving the element routine the residual and tangent are modified for the Dirichlet condition.  The Dirichlet equations are not assembled in the global problem.

ParFEAP can do either Option 1 or Option 2:

[Prof. S. Govindjee]

To clarify, when using parFEAP the two options are (a) solve a reduced size system and (b) solve a full size system.  When (b) is chosen, the so-called boundary conditions in method, then FEAP inserts a 1 on the diagonal (zeroed out rows and columns) for the dirichlet dofs, with the boundary value inserted on the RHS.

[festein]

First, thank you for your answers, but they confuse me a bit. The parfeap implementation is not important to me, since I only solve small problems.
You said "After leaving the element routine the residual and tangent are modified for the Dirichlet condition." What do you mean by that? Reducinig the system or changing the values of the tangent? I don't know how this would be done, since the element routine defines the properties of the tangent.
In my understanding the procedure looks as follows:
- Initialize u_glob (without dirichlet boundary conditions)
- Create the element stiffness matrices k_elem and element residuals r_elem
- Assemble K_glob and R_glob based on k_elem and r_elem; disregard the rows/columns where a Dirichlet boundary condition is applied
- Modify R_glob based on the dirichlet values
- solve the system
Am i correct?
Of course if I got a simulation involving multiple steps (for example J2Plasticity) the initialization of u_glob would be the one from the last time step.

[Prof. S. Govindjee]

The relevant routines are in the program/ folder.  The flow is roughly

formfe -> pform -> plocal -> elmlib -> modify -> dasble

formfe initiates the forming of the arrays.  pform loops over the elements, calling plocal which localizes global values of X and U into local arrays XL and UL (for all dofs on the element), then elemlib is called to get the element residual, p_elem, (Btranspose sigma + body forces) and stiffness, s_elem, then modify is called to adjust the values in p_elem by the columns of s_elem times and driven boundary increments, finally the rows/columns of p_elem and s_elem are assembled in to the global stiffness, during the assembly only the active dofs are assembled; dofs with driven values are not assembled.  Note that at this stage the local element residual has already been modified; i.e. the modification for the driven dofs is performed at the local level, before assembly.

[JStorm]

Which means, that the assembly of the element matrices is towards the reduced global system of equations.
The original system of questions is never constructed.

[Prof. S. Govindjee]

Correct the full system is not constructed.  The closest one gets to that is when one wants to compute reactions.

[festein]

Thank you, that was exactly what I needed to know. Maybe it would be useful to put that into the theory manual as well?