fem

package module
v0.0.0-...-36c816c Latest Latest
Warning

This package is not in the latest version of its module.

 Go to latest Published: Jul 9, 2023 License: MIT Imports: 8 Imported by: 0

README

go-fem

Extremely clean API toolkit for finite element analysis/methods.

Roadmap
  • 3D Isoparametric element assembly.
  • Better sparse assembler As 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

Examples

Constants

View Source 
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 

type Constituter interface {
	Constitutive() (mat.Matrix, error)
}

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
)

func (DofsFlag) Count 

func (d DofsFlag) Count() int

Count returns the number of dofs set in d.

func (DofsFlag) Has 

func (d DofsFlag) Has(q DofsFlag) bool

Has returns true if d has all of q's dofs set. It returns false if q has a dof set that is not set in d.

func (DofsFlag) String 

func (d DofsFlag) String() (s string)

String returns a human readable representation of which dofs are set in d.

type Element 

type Element interface {
	// LenNodes returns the number of nodes in the element.
	LenNodes() int
	// Dofs returns the degrees of freedom corresponding to each of the element's nodes.
	Dofs() DofsFlag
}

type Element3 

type Element3 interface {
	Element
	CopyK(dst *mat.Dense, elementNodes []r3.Vec) error
	SetConstitutive(c Constituter) error
}

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 NewFixity 

func NewFixity(modelDofs DofsFlag, numNodes int) Fixity

NewFixity returns a new Fixity with the given modelDofs and number of nodes.

func (Fixity) Fix 

func (f Fixity) Fix(nodeIdx int, fixed DofsFlag)

Fix sets the fixed dofs for the given node. It ignores dofs that are not in the model.

func (Fixity) FixedDofs 

func (f Fixity) FixedDofs() []int

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.

func (Fixity) Free 

func (f Fixity) Free(nodeIdx int, freed DofsFlag)

Free sets the free dofs for the given node. It ignores dofs that are not in the model.

func (Fixity) FreeDofs 

func (f Fixity) FreeDofs() []int

FreeDofs returns the free 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. 

elemType := elements.Tetra4{}
steel := solids.Isotropic{E: 200e9, Poisson: 0.3}
nod := []r3.Vec{
	{X: 0, Y: 0, Z: 0},
	{X: 1}, // is (1,0,0) point
	{Y: 1},
	{Z: 1},
}
elems := [][]int{
	{0, 1, 2, 3}, // one single element.
}
ga := fem.NewGeneralAssembler(nod, elemType.Dofs())
err := ga.AddIsoparametric(elemType, steel.Solid3D(), len(elems), func(i int) ([]int, r3.Vec, r3.Vec) {
	return elems[i], r3.Vec{}, r3.Vec{}
})
if err != nil {
	panic(err)
}
fmt.Printf("K=\n%.5g", lap.Formatted(ga.Ksolid()))

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) 
// 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⎦
Example (AxisymmetricThermal) 
// We assemble a single 2D axisymmetric element.
materialTherm := solids.IsotropicConductivity{K: 45} // 45 W/mK for steel.
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{NodeDofs: fem.DofPosX}
ga := fem.NewGeneralAssembler(nodes, elemType.Dofs())
err := ga.AddIsoparametric(elemType, materialTherm.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=
⎡  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) 
// We assemble a single Quad4 element with a unitary isotropic material.
material := solids.Isotropic{E: 1, Poisson: 0.3}
nodes := []r3.Vec{
	{X: 0, Y: 0},
	{X: 1, Y: 0},
	{X: 1, Y: 1},
	{X: 0, Y: 1},
}
elemType := elements.Quad4{}
elems := [][]int{
	{0, 1, 2, 3},
}

ga := fem.NewGeneralAssembler(nodes, elemType.Dofs())
err := ga.AddIsoparametric(elemType, material.PlaneStess(), len(elems), func(i int) (elem []int, xC, yC r3.Vec) {
	return elems[i], r3.Vec{}, r3.Vec{}
})
if err != nil {
	panic(err)
}
fmt.Printf("K=\n%.3f", lap.Formatted(ga.Ksolid()))

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 IsoConstituter interface {
	Constituter
	SetStrainDisplacementMatrix(dstB, elemNodes, dN *mat.Dense, N *mat.VecDense) (scale float64)
}

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.

Source Files

View all Source files

Directories

Path Synopsis
constitution
solids
package solids provides constitutive models for solid mechanics analysis.
 package solids provides constitutive models for solid mechanics analysis.
examples
axisymmetric
axisymmetric_thermal
bimetallic
ruc
internal
expmat

Jump to

Keyboard shortcuts

? : This menu
/ : Search site
f or F : Jump to
y or Y : Canonical URL
go.dev uses cookies from Google to deliver and enhance the quality of its services and to analyze traffic. Learn more.