#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------
#
# _____ _ _ _ _____ _____
# | ____|__| | ___| |_ _____(_)___ ___| ___| ____|
# | _| / _` |/ _ \ \ \ /\ / / _ \ / __/ __| |_ | _|
# | |__| (_| | __/ |\ V V / __/ \__ \__ \ _| | |___
# |_____\__,_|\___|_| \_/\_/ \___|_|___/___/_| |_____|
#
#
# Unit of Strength of Materials and Structural Analysis
# University of Innsbruck,
# 2017 - today
#
# Matthias Neuner matthias.neuner@uibk.ac.at
#
# This file is part of EdelweissFE.
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License as published by the Free Software Foundation; either
# version 2.1 of the License, or (at your option) any later version.
#
# The full text of the license can be found in the file LICENSE.md at
# the top level directory of EdelweissFE.
# ---------------------------------------------------------------------
# Created on Fri Aug 5 08:35:06 2022
# @author: Matthias Neuner
"""
Implementing your own finite elements can be done easily by subclassing from
the abstract base class :class:`~BaseElement`.
"""
from abc import abstractmethod
import numpy as np
from edelweissfe.nodecouplingentity.base.nodecouplingentity import (
BaseNodeCouplingEntity,
)
from edelweissfe.numerics.vijentitybase import VIJEntityBase
from edelweissfe.points.node import Node
[docs]
class BaseElement(BaseNodeCouplingEntity, VIJEntityBase):
@property
@abstractmethod
def elNumber(self) -> int:
"""The unique number of this element"""
@property
@abstractmethod
def elType(self) -> str:
"""The type of this element."""
@property
@abstractmethod
def hasMaterial(self) -> str:
"""Flag to check if a material was assigned to this element."""
@abstractmethod
def __init__(self, elementType: str, elNumber: int):
"""Finite elements in EdelweissFE should be derived from this
base class in order to follow the general interface.
Every method that computes a residual contribution ``P`` with its tangent ``K = dP/dU``
follows one rule (see :ref:`residual-and-tangent-convention`): :meth:`computeKernels` adds
the internal forces, the load methods :meth:`computeDistributedLoad` and
:meth:`computeBodyForce` add the *negative* of the external load. EdelweissFE expects the
layout of ``P`` (and of ``K``) to be of the form
.. code-block:: console
[ node 1 - dofs field 1,
node 1 - dofs field 2,
node 1 - ... ,
node 1 - dofs field n,
node 2 - dofs field 1,
... ,
node N - dofs field n].
Parameters
----------
elementType
A string identifying the requested element formulation.
elNumber
A unique integer label used for all kinds of purposes.
"""
[docs]
@abstractmethod
def setNodes(self, nodes: list[Node]):
"""Assign the nodes to the element.
Parameters
----------
nodes
A list of nodes.
"""
[docs]
@abstractmethod
def setProperties(self, elementProperties: np.ndarray):
"""Assign a set of properties to the element.
Parameters
----------
elementProperties
A numpy array containing the element properties.
"""
[docs]
def assignProperty(self, propertyName: str, properties: np.ndarray):
"""Assign a property of the element by name.
Parameters
----------
propertyName
The name of the property.
properties
A numpy array containing the property values.
"""
raise NotImplementedError("This element type does not support named property assignment.")
[docs]
def getPropertyNames(self) -> list[str]:
"""Get the names of all the valid properties of the element.
Returns
-------
list[str]
A list of valid property names.
"""
return []
@property
def propertyNames(self) -> list[str]:
"""Get the names of all the valid properties of the element.
Returns
-------
list[str]
A list of valid property names.
"""
return self.getPropertyNames()
[docs]
@abstractmethod
def initializeElement(
self,
):
"""Initalize the element to be ready for computing."""
[docs]
@abstractmethod
def setMaterial(self, materialName: str, materialProperties: np.ndarray):
"""Assign a material and material properties.
Parameters
----------
materialName
The name of the requested material.
materialProperties
The numpy array containing the material properties.
"""
[docs]
def updateMaterialProperty(self, index: int, value: float):
"""Change one entry of the property vector of the assigned material, during an analysis.
Unlike :meth:`setMaterial`, this keeps the state of the element: the state variables (e.g.
stress, plastic strain, damage) are neither reset nor re-initialized, and the arrays holding
them stay the same objects, so persistent result views (field outputs) remain valid. All
other properties, including any assigned per element (e.g. ``materialParameterFromField``),
are kept.
Elements without a material do not support this.
Parameters
----------
index
The index of the property in the material's property vector.
value
The new value of the property.
"""
raise NotImplementedError(f"{type(self).__name__} does not support changing the properties of its material.")
[docs]
@abstractmethod
def setInitialCondition(self, stateType: str, values: np.ndarray):
"""Assign initial conditions.
Parameters
----------
stateType
The type of initial state.
values
The numpy array describing the initial state.
"""
[docs]
@abstractmethod
def computeDistributedLoad(
self,
loadType: str,
P: np.ndarray,
K: np.ndarray,
faceID: int,
load: np.ndarray,
U: np.ndarray,
time: float,
dT: float,
):
"""Add the residual contribution of a surface load, the negative of the external load
:math:`-\\boldsymbol{f}_\\mathrm{ext}`, and its tangent, which is nonzero only for a load that
follows the deformation (e.g. a pressure on the current surface).
Parameters
----------
loadType
The type of load.
P
The residual contribution to augment, :math:`-\\boldsymbol{f}_\\mathrm{ext}`.
K
The tangent ``dP/dU`` to augment.
faceID
The number of the elements face this load acts on.
load
The magnitude (or vector) describing the load.
U
The current solution vector.
time
The current time.
dTime
The time increment.
"""
[docs]
@abstractmethod
def computeKernels(
self,
K: np.ndarray,
P: np.ndarray,
U: np.ndarray,
dU: np.ndarray,
time: float,
dT: float,
):
"""Add the internal forces to the residual, and their tangent, for given time, field, and field
increment.
Parameters
----------
K
The tangent ``dP/dU`` to augment, the stiffness matrix.
P
The residual contribution to augment, the internal forces :math:`\\boldsymbol{f}_\\mathrm{int}`.
U
The current solution vector.
dU
The current solution vector increment.
time
The current time.
dTime
The time increment.
"""
[docs]
@abstractmethod
def computeKernelsExplicit(
self,
P: np.ndarray,
U: np.ndarray,
dU: np.ndarray,
time: float,
dT: float,
):
"""Add the internal forces to the residual for given time, field, and field increment, without a
tangent.
Parameters
----------
P
The residual contribution to augment, the internal forces :math:`\\boldsymbol{f}_\\mathrm{int}`.
U
The current solution vector.
dU
The current solution vector increment.
time
The current time.
dTime
The time increment.
"""
[docs]
@abstractmethod
def computeLumpedInertia(
self,
M: np.ndarray,
):
"""Evaluate the diagonal of the lumped inertia of the element, over every field it carries.
The coefficient of each field's SECOND time derivative: mass on the displacement block,
and a micro-inertia on a non-local block that has been given one. A parabolic
gradient-enhanced field reports none, which is what tells the solver it is first order in
time; see computeLumpedDamping() for what integrates it then.
Parameters
----------
M
The diagonal of the lumped inertia to be defined.
"""
[docs]
def computeConsistentInertia(
self,
M: np.ndarray,
):
"""Evaluate the consistent (full) inertia matrix of the element, over every field it
carries.
The counterpart of :meth:`computeLumpedInertia` for an implicit dynamic solver, which
assembles a sparse global mass matrix rather than dividing by a diagonal: the coefficient
matrix of each field's SECOND time derivative, :math:`\\int_\\Omega \\rho N^T N \\, dV` on
the displacement block. Written into ``M`` in the same layout the element uses for its
stiffness in ``computeKernels``, so that a slice of the VIJ system matrix can be handed over
directly and the assembled mass shares the stiffness' sparsity pattern.
Not abstract: an element formulation without one is refused at assembly time by the solver
that needs it, rather than every element being forced to implement a matrix only the
implicit dynamic solver reads.
Parameters
----------
M
The consistent inertia to be defined, in the element's stiffness layout.
Raises
------
NotImplementedError
If the element does not provide a consistent inertia.
"""
raise NotImplementedError(
"{:} does not provide a consistent inertia matrix (computeConsistentInertia).".format(type(self).__name__)
)
[docs]
def computeLumpedDamping(
self,
C: np.ndarray,
):
"""Evaluate the diagonal of the lumped damping of the element, over every field it carries.
The coefficient of each field's FIRST time derivative: zero on the displacement block,
and the non-local viscosity on a non-local block -- always, whether or not that field also
carries a micro-inertia. A first-order field is integrated by this term alone; a
second-order one is damped by it. Carrying none is the ordinary case, so the default
leaves the buffer alone.
Parameters
----------
C
The diagonal of the lumped damping, supplied zero-initialised and written in place.
Entries left untouched are read as "no damping here".
"""
@property
def initialVelocity(self) -> np.ndarray:
"""The element's initial velocity, in DOF order, applied once as an
initial condition at the start of an explicit simulation.
Elements do not take part in the per-step momentum remap (that is a
particle/cell concern); a mass-bearing element such as
:class:`~edelweissfe.elements.pointmass.PointMass` instead declares its
initial velocity here, which the explicit solver seeds into the grid
velocity vector once. The default is rest.
Returns
-------
np.ndarray
The initial velocity of the element's DOFs (size ``nDof``).
"""
return np.zeros(self.nDof)
@property
def hasKernels(self) -> bool:
"""Whether this entity contributes to the internal force and carries a state of its own.
Not every entity in the element container is a finite element. A contact facet is kept
there for its nodes and its geometry alone: it has no material and no state, and its
kernels, its internal energy and its state acceptance are all no-ops. A solver that
touches every element on every increment may leave those out instead of calling into
them for nothing.
Returns
-------
bool
True for a finite element proper, False for an entity that carries geometry only.
"""
return True
[docs]
@abstractmethod
def computeCriticalTimeStepForExplicitDynamics(
self,
Q: np.ndarray,
) -> float:
"""Evaluate the critical time step for explicit dynamics.
Parameters
----------
Q
The current solution vector, which might be needed to compute the critical time step due to nonlinearities.
Returns
-------
float
The critical time step.
"""
[docs]
@abstractmethod
def computeInternalEnergy(self) -> float:
"""Evaluate the internal energy of the element.
Returns
-------
float
The internal energy.
"""
[docs]
@abstractmethod
def computeBodyForce(
self,
P: np.ndarray,
K: np.ndarray,
load: np.ndarray,
U: np.ndarray,
time: float,
dTime: float,
):
"""Add the residual contribution of a body force, the negative of the external load
:math:`-\\boldsymbol{f}_\\mathrm{ext}`, and its tangent, which is zero for a dead load.
Parameters
----------
P
The residual contribution to augment, :math:`-\\boldsymbol{f}_\\mathrm{ext}`.
K
The tangent ``dP/dU`` to augment.
load
The magnitude (or vector) describing the load.
U
The current solution vector.
time
The current time.
dTime
The time increment.
"""
[docs]
@abstractmethod
def acceptLastState(
self,
):
"""Accept the computed state (in nonlinear iteration schemes)."""
[docs]
@abstractmethod
def resetToLastValidState(
self,
):
"""Rest to the last valid state."""
[docs]
def getStateVars(self) -> np.ndarray:
"""Return a copy of the element's flat (converged) quadrature-point state-variable buffer.
The buffer is laid out as ``nQuadraturePoints`` contiguous per-point blocks. Used by adaptive
mesh refinement to hand history down to child elements. Override in concrete elements that
carry a state buffer.
Returns
-------
np.ndarray
A copy of the converged state-variable buffer.
"""
raise NotImplementedError("This element does not expose a state-variable buffer.")
[docs]
def setStateVars(self, values: np.ndarray):
"""Overwrite the element's (converged and trial) state-variable buffer.
Parameters
----------
values
The new state-variable buffer, same size as :meth:`getStateVars`.
"""
raise NotImplementedError("This element does not expose a state-variable buffer.")
[docs]
def getStateVarSlice(self, name: str) -> tuple:
"""Locate a named state variable within a single per-quadrature-point state block.
Enables per-state-variable routing of adaptive-refinement state transfer (see
:class:`~edelweissfe.adaptivity.statetransfer.perstatevar.PerStateVarStateTransfer`).
Parameters
----------
name
The state-variable name (as understood by the underlying material / element).
Returns
-------
tuple
``(offset, size)`` of the named variable within one per-quadrature-point block.
"""
raise NotImplementedError("This element does not expose named state-variable slices.")
[docs]
@abstractmethod
def getResultArray(self, result: str, quadraturePoint: int, getPersistentView: bool = True) -> np.ndarray:
"""Get the array of a result, possibly as a persistent view which is continiously
updated by the element.
Parameters
----------
result
The name of the result.
quadraturePoint
The number of the quadrature point.
getPersistentView
If true, the returned array should be continiously updated by the element.
Returns
-------
np.ndarray
The result.
"""
[docs]
@abstractmethod
def getCoordinatesAtCenter(self) -> np.ndarray:
"""Compute the underlying MarmotElement centroid coordinates.
Returns
-------
np.ndarray
The element's central coordinates.
"""
[docs]
@abstractmethod
def getCoordinatesAtQuadraturePoints(self) -> np.ndarray:
"""Compute the underlying MarmotElement qp coordinates.
Returns
-------
np.ndarray
The element's qp coordinates.
"""
[docs]
@abstractmethod
def getNumberOfQuadraturePoints(self) -> int:
"""Compute the underlying MarmotElement qp coordinates.
Returns
-------
np.ndarray
The element's qp coordinates.
"""