README
¶
go-fem
Extremely clean API toolkit for finite element analysis/methods.
Roadmap
3D Isoparametric element assembly.Better sparse assemblerAs fast as high performance libraries.2D Isoparametric element assembly.Arbitrary Isoparametric element assembly.- Stress extraction from displacements.
- Shell and plate element assembly.
- Better define Element3 API.
- Add constraints assembly.
- Lagrange Multipliers or Penalty method.
- Stiffness matrix utility functions for better conditioning and troubleshooting.
Axisymmetric assembly example
package main
import (
"fmt"
"log"
"github.com/soypat/go-fem"
"github.com/soypat/go-fem/constitution/solids"
"github.com/soypat/go-fem/elements"
"github.com/soypat/lap"
"gonum.org/v1/gonum/spatial/r3"
)
func main() {
// We assemble a single 2D axisymmetric element.
// With a young modulus of 1000 and a poisson ratio of 0.33.
// The units are undefined. If using SI, E would be 1kPa and node positions
// would be in meters.
material := solids.Isotropic{E: 1000, Poisson: 0.33}
nodes := []r3.Vec{
{X: 20, Y: 0},
{X: 30, Y: 0},
{X: 30, Y: 1},
{X: 20, Y: 1},
}
elems := [][4]int{
{0, 1, 2, 3},
}
elemType := elements.Quad4{}
ga := fem.NewGeneralAssembler(nodes, elemType.Dofs())
err := ga.AddIsoparametric(elemType, material.Axisymmetric(), 1, func(i int) (elem []int, xC, yC r3.Vec) {
return elems[i][:], r3.Vec{}, r3.Vec{}
})
if err != nil {
log.Fatal(err)
}
fmt.Printf("K=\n%.3g\n", lap.Formatted(ga.Ksolid()))
}
Output:
K=
⎡ 2.93e+04 5.23e+03 1.45e+04 2.06e+03 -1.63e+04 -6.45e+03 -2.77e+04 -848⎤
⎢ 5.23e+03 1.11e+05 -2.36e+03 6.14e+04 -7.37e+03 -6.19e+04 848 -1.11e+05⎥
⎢ 1.45e+04 -2.36e+03 3.6e+04 -8.59e+03 -3.37e+04 3.58e+03 -1.63e+04 7.37e+03⎥
⎢ 2.06e+03 6.14e+04 -8.59e+03 1.36e+05 -3.58e+03 -1.36e+05 6.45e+03 -6.19e+04⎥
⎢-1.63e+04 -7.37e+03 -3.37e+04 -3.58e+03 3.6e+04 8.59e+03 1.45e+04 2.36e+03⎥
⎢-6.45e+03 -6.19e+04 3.58e+03 -1.36e+05 8.59e+03 1.36e+05 -2.06e+03 6.14e+04⎥
⎢-2.77e+04 848 -1.63e+04 6.45e+03 1.45e+04 -2.06e+03 2.93e+04 -5.23e+03⎥
⎣ -848 -1.11e+05 7.37e+03 -6.19e+04 2.36e+03 6.14e+04 -5.23e+03 1.11e+05⎦
Benchmarks
Takes ~33s to assemble 680k dofs for 1.3 million 4 node 3D tetrahedrons on a low end AMD Zen2 CPU:
$ go test -bench=. -benchmem .
goos: linux
goarch: amd64
pkg: github.com/soypat/go-fem
cpu: AMD Ryzen 5 3400G with Radeon Vega Graphics
BenchmarkTetra4Assembly/105_dofs,_48_elems-8 1218 932708 ns/op 155646 B/op 744 allocs/op
BenchmarkTetra4Assembly/3723_dofs,_5376_elems-8 10 109067373 ns/op 10928022 B/op 74500 allocs/op
BenchmarkTetra4Assembly/27027_dofs,_46080_elems-8 1 1023632489 ns/op 88685256 B/op 634443 allocs/op
BenchmarkTetra4Assembly/206115_dofs,_380928_elems-8 1 8849933029 ns/op 717490648 B/op 5229747 allocs/op
BenchmarkTetra4Assembly/684723_dofs,_1299456_elems-8 1 32997698741 ns/op 2742313160 B/op 17888792 allocs/op
Documentation
¶
Index ¶
- Constants
- type Constituter
- type DofsFlag
- type Element
- type Element3
- type ElementInfo
- type Fixity
- type GeneralAssembler
- func (ga *GeneralAssembler) AddElement3(elemT Element3, c Constituter, Nelem int, ...) error
- func (ga *GeneralAssembler) AddIsoparametric(elemT Isoparametric, c IsoConstituter, Nelem int, ...) error
- func (ga *GeneralAssembler) DofMapping(e Element) ([]int, error)
- func (ga *GeneralAssembler) ElementInfo(elemT Element, element []int) (einfo ElementInfo, _ error)
- func (ga *GeneralAssembler) ForEachElement(elemT Element, spatialDimsPerNode, Nelem int, getElement func(i int) []int, ...) error
- func (ga *GeneralAssembler) IsoparametricStrains(displacements lap.Vector, elemT Isoparametric, c IsoConstituter, Nelem int, ...) error
- func (ga *GeneralAssembler) Ksolid() *lap.Sparse
- func (ga *GeneralAssembler) TotalDofs() int
- type IsoConstituter
- type Isoparametric
Examples ¶
Constants ¶
const ( DofPos = DofPosX | DofPosY | DofPosZ DofRot = DofRotX | DofRotY | DofRotZ Dof6 = DofPos | DofRot )
Common degree of freedom flag definitions.
Variables ¶
This section is empty.
Functions ¶
This section is empty.
Types ¶
type Constituter ¶
Constituter represents the homogenous properties of a medium that can then be used to model solids or other continuous field problems. For solids it returns the unmodified constitutive tensor (Generalized Hookes law). The shear modulii should not be halved.
type DofsFlag ¶
type DofsFlag uint16
DofsFlag holds bitwise information of degrees of freedom.
const ( // DofPosX corresponds to X position degree of freedom (0). DofPosX DofsFlag = 1 << iota // DofPosY corresponds to Y position degree of freedom (1). DofPosY // DofPosZ corresponds to Z position degree of freedom (2). DofPosZ // DofRotX corresponds to rotational degree of freedom about X axis (3). DofRotX // DofRotY corresponds to rotational degree of freedom about Y axis (4). DofRotY // DofRotZ corresponds to rotational degree of freedom about Z axis (5). DofRotZ )
type ElementInfo ¶
type ElementInfo struct {
// contains filtered or unexported fields
}
func (ElementInfo) Dofs ¶
func (einfo ElementInfo) Dofs() [][]int
Dofs returns the element's dofs:
- The first index "n" is the node index.
- The second index "d" is corresponds to the dth dof of the nth node.
func (ElementInfo) Nodes ¶
func (einfo ElementInfo) Nodes() []r3.Vec
Nodes returns the element's nodes in order of the element int slice used to create the ElementInfo.
type Fixity ¶
type Fixity struct {
// contains filtered or unexported fields
}
Fixity holds information for the fixed degrees of freedom per node. It assumes constant number of dofs per node.
func (Fixity) Fix ¶
Fix sets the fixed dofs for the given node. It ignores dofs that are not in the model.
func (Fixity) FixedDofs ¶
FixedDofs returns the fixed dof indices. This should be used to set the free dofs in the global stiffness matrix using lap.Slice or similar slicing functions.
type GeneralAssembler ¶
type GeneralAssembler struct {
// contains filtered or unexported fields
}
GeneralAssembler is a general purpose stiffness matrix assembler.
Example ¶
Calculates the stiffness matrix of a single isoparametric tetrahedron element.
Output: K= ⎡ 4.2308e+11 1.9231e+11 1.9231e+11 -2.6923e+11 -7.6923e+10 -7.6923e+10 -7.6923e+10 -1.1538e+11 0 -7.6923e+10 0 -1.1538e+11⎤ ⎢ 1.9231e+11 4.2308e+11 1.9231e+11 -1.1538e+11 -7.6923e+10 0 -7.6923e+10 -2.6923e+11 -7.6923e+10 0 -7.6923e+10 -1.1538e+11⎥ ⎢ 1.9231e+11 1.9231e+11 4.2308e+11 -1.1538e+11 0 -7.6923e+10 0 -1.1538e+11 -7.6923e+10 -7.6923e+10 -7.6923e+10 -2.6923e+11⎥ ⎢-2.6923e+11 -1.1538e+11 -1.1538e+11 2.6923e+11 0 0 0 1.1538e+11 0 0 0 1.1538e+11⎥ ⎢-7.6923e+10 -7.6923e+10 0 0 7.6923e+10 0 7.6923e+10 0 0 0 0 0⎥ ⎢-7.6923e+10 0 -7.6923e+10 0 0 7.6923e+10 0 0 0 7.6923e+10 0 0⎥ ⎢-7.6923e+10 -7.6923e+10 0 0 7.6923e+10 0 7.6923e+10 0 0 0 0 0⎥ ⎢-1.1538e+11 -2.6923e+11 -1.1538e+11 1.1538e+11 0 0 0 2.6923e+11 0 0 0 1.1538e+11⎥ ⎢ 0 -7.6923e+10 -7.6923e+10 0 0 0 0 0 7.6923e+10 0 7.6923e+10 0⎥ ⎢-7.6923e+10 0 -7.6923e+10 0 0 7.6923e+10 0 0 0 7.6923e+10 0 0⎥ ⎢ 0 -7.6923e+10 -7.6923e+10 0 0 0 0 0 7.6923e+10 0 7.6923e+10 0⎥ ⎣-1.1538e+11 -1.1538e+11 -2.6923e+11 1.1538e+11 0 0 0 1.1538e+11 0 0 0 2.6923e+11⎦
func NewGeneralAssembler ¶
func NewGeneralAssembler(nodes []r3.Vec, modelDofs DofsFlag) *GeneralAssembler
NewSymAssembler initializes a GeneralAssembler ready for use.
func (*GeneralAssembler) AddElement3 ¶
func (ga *GeneralAssembler) AddElement3(elemT Element3, c Constituter, Nelem int, getElement func(i int) (e []int, x, y r3.Vec)) error
func (*GeneralAssembler) AddIsoparametric ¶
func (ga *GeneralAssembler) AddIsoparametric(elemT Isoparametric, c IsoConstituter, Nelem int, getElement func(i int) (elem []int, xC, yC r3.Vec)) error
AddIsoparametric adds isoparametric elements to the model's solid stiffness matrix. It calls getElement to get the element nodes and coordinates Nelem times, each with an incrementing element index i.
Arbitrary orientation of solid properties for each isoparametric element is not yet implemented.
Example (Axisymmetric) ¶
Output: K= ⎡ 2.93e+04 5.23e+03 1.45e+04 2.06e+03 -1.63e+04 -6.45e+03 -2.77e+04 -848⎤ ⎢ 5.23e+03 1.11e+05 -2.36e+03 6.14e+04 -7.37e+03 -6.19e+04 848 -1.11e+05⎥ ⎢ 1.45e+04 -2.36e+03 3.6e+04 -8.59e+03 -3.37e+04 3.58e+03 -1.63e+04 7.37e+03⎥ ⎢ 2.06e+03 6.14e+04 -8.59e+03 1.36e+05 -3.58e+03 -1.36e+05 6.45e+03 -6.19e+04⎥ ⎢-1.63e+04 -7.37e+03 -3.37e+04 -3.58e+03 3.6e+04 8.59e+03 1.45e+04 2.36e+03⎥ ⎢-6.45e+03 -6.19e+04 3.58e+03 -1.36e+05 8.59e+03 1.36e+05 -2.06e+03 6.14e+04⎥ ⎢-2.77e+04 848 -1.63e+04 6.45e+03 1.45e+04 -2.06e+03 2.93e+04 -5.23e+03⎥ ⎣ -848 -1.11e+05 7.37e+03 -6.19e+04 2.36e+03 6.14e+04 -5.23e+03 1.11e+05⎦
Example (AxisymmetricThermal) ¶
Output: K= ⎡ 3.4e+03 1.84e+03 -1.89e+03 -3.36e+03⎤ ⎢ 1.84e+03 4.18e+03 -4.1e+03 -1.89e+03⎥ ⎢-1.89e+03 -4.1e+03 4.18e+03 1.84e+03⎥ ⎣-3.36e+03 -1.89e+03 1.84e+03 3.4e+03⎦
Example (Quad4PlaneStress) ¶
Output: K= ⎡ 0.495 0.179 -0.302 -0.014 -0.247 -0.179 0.055 0.014⎤ ⎢ 0.179 0.495 0.014 0.055 -0.179 -0.247 -0.014 -0.302⎥ ⎢-0.302 0.014 0.495 -0.179 0.055 -0.014 -0.247 0.179⎥ ⎢-0.014 0.055 -0.179 0.495 0.014 -0.302 0.179 -0.247⎥ ⎢-0.247 -0.179 0.055 0.014 0.495 0.179 -0.302 -0.014⎥ ⎢-0.179 -0.247 -0.014 -0.302 0.179 0.495 0.014 0.055⎥ ⎢ 0.055 -0.014 -0.247 0.179 -0.302 0.014 0.495 -0.179⎥ ⎣ 0.014 -0.302 0.179 -0.247 -0.014 0.055 -0.179 0.495⎦
func (*GeneralAssembler) DofMapping ¶
func (ga *GeneralAssembler) DofMapping(e Element) ([]int, error)
DofMapping creates element to model dof mapping. If model does not contain all argument elements dofs then an error is returned. The length of the slice returned is equal to the amount of dofs per element node.
func (*GeneralAssembler) ElementInfo ¶
func (ga *GeneralAssembler) ElementInfo(elemT Element, element []int) (einfo ElementInfo, _ error)
ElementInfo returns the element's nodes and dofs in order of the element int slice. The Dofs corresponde
func (*GeneralAssembler) ForEachElement ¶
func (ga *GeneralAssembler) ForEachElement(elemT Element, spatialDimsPerNode, Nelem int, getElement func(i int) []int, elemCallback elementDofCallback) error
IntegrateIsoparametric is a general purpose for-each isoparametric element iterator. This function iterates over Nelem elements of type elemT, calling getElement to get the element nodes indices. These are then used to get the element's nodes and dofs which are passed in the elemCallback function.
func (*GeneralAssembler) IsoparametricStrains ¶
func (ga *GeneralAssembler) IsoparametricStrains(displacements lap.Vector, elemT Isoparametric, c IsoConstituter, Nelem int, getElement func(i int) (elem []int, xC, yC r3.Vec), strainCallback func(iele int, strains []float64)) error
IsoparametricStrains calculates the strains at the integration points of an isoparametric element.
func (*GeneralAssembler) Ksolid ¶
func (ga *GeneralAssembler) Ksolid() *lap.Sparse
Ksolid returns the stiffness matrix of the solid.
func (*GeneralAssembler) TotalDofs ¶
func (ga *GeneralAssembler) TotalDofs() int
TotalDofs returns the total number of dofs in the model.
type IsoConstituter ¶
type Isoparametric ¶
type Isoparametric interface { Element // BasisNodes returns the positions of the nodes relative to the origin // of the element. Returned slice is of length LenNodes. IsoparametricNodes() []r3.Vec // Quadrature returns the desired quadrature integration positions and // respective weights. Quadrature() (positions []r3.Vec, weights []float64) // Basis returns the form functions of the element evaluated at a position // relative to the origin of the element. Returned slice is of length LenNodes. Basis(r3.Vec) []float64 // BasisDiff returns the differentiated form functions of the element // evaluated at a position relative to the origin of the element. // The result first contains the form functions differentiated // with respect to X, then Y and finally Z. // [ dN/dx ] // N = | dN/dy | (row major form) // [ dN/dz ] // Suggested way of initializing matrix: // dN := mat.NewDense(3, e.LenNodes(), e.BasisDiff(v)) // Returned slice is of length LenNodes*3. BasisDiff(r3.Vec) []float64 }
Isoparametric is a 3D/2D/1D element that employs the isoparametric formulation.