fem

type Beam struct {

	// basic data
	Cid int         // cell/element id
	X   [][]float64 // matrix of nodal coordinates [ndim][nnode]
	Nu  int         // total number of unknowns == 2 * nsn

	// parameters
	E   float64 // Young's modulus
	A   float64 // cross-sectional area
	Izz float64 // Inertia zz

	// variables for dynamics
	Rho  float64  // density of solids
	Gfcn fun.Func // gravity function

	// vectors and matrices
	T   [][]float64 // global-to-local transformation matrix [nnode*ndim][nnode*ndim]
	Kl  [][]float64 // local K matrix
	K   [][]float64 // global K matrix
	Ml  [][]float64 // local M matrices
	M   [][]float64 // global M matrices
	Rus []float64   // residual: Rus = fi - fx

	// problem variables
	Umap []int    // assembly map (location array/element equations)
	Hasq bool     // has distributed loads
	QnL  fun.Func // distributed normal load functions: left
	QnR  fun.Func // distributed normal load functions: right
	Qt   fun.Func // distributed tangential load
	// contains filtered or unexported fields
}

type Decoder interface {
	Decode(e interface{}) error
}

type Dof struct {
	Key string // primary variable key. e.g. "ux"
	Eq  int    // equation number
}

type Domain struct {

	// init: region, mesh, linear solver
	Reg    *inp.Region // region data
	Msh    *inp.Mesh   // mesh data
	LinSol la.LinSol   // linear solver

	// stage: auxiliary maps for setting boundary conditions
	FaceConds map[int][]*FaceCond // maps cell id to its face boundary conditions

	// stage: nodes (active) and elements (active AND in this processor)
	Nodes  []*Node // active nodes (for each stage)
	Elems  []Elem  // [procNcells] only active elements in this processor (for each stage)
	MyCids []int   // [procNcells] the ids of cells in this processor

	// stage: auxiliary maps for dofs and equation types
	F2Y      map[string]string // converts f-keys to y-keys; e.g.: "ux" => "fx"
	YandC    map[string]bool   // y and constraints keys; e.g. "ux", "pl", "H", "incsup", "rigid"
	Dof2Tnum map[string]int    // {t1,t2}-types: dof => t_number; e.g. "ux" => 2, "pl" => 1

	// stage: auxiliary maps for nodes and elements
	Vid2node   []*Node // [nverts] VertexId => index in Nodes. Inactive vertices are 'nil'
	Cid2elem   []Elem  // [ncells] CellId => index in Elems. Cells in other processors or inactive are 'nil'
	Cid2active []bool  // [ncells] CellId => whether cell is active or not in ANY processor

	// stage: subsets of elements
	ElemIntvars []ElemIntvars   // elements with internal vars in this processor
	ElemConnect []ElemConnector // connector elements in this processor

	// stage: coefficients and prescribed forces
	EssenBcs EssentialBcs // constraints (Lagrange multipliers)
	PtNatBcs PtNaturalBcs // point loads such as prescribed forces at nodes

	// stage: t1 and t2 variables
	T1eqs []int // first t-derivative variables; e.g.:  dp/dt vars (subset of ykeys)
	T2eqs []int // second t-derivative variables; e.g.: d²u/dt² vars (subset of ykeys)

	// stage: dimensions
	NnzKb int // number of nonzeros in Kb matrix
	Ny    int // total number of dofs, except λ
	Nlam  int // total number of Lagrange multipliers
	NnzA  int // number of nonzeros in A (constraints) matrix
	Nyb   int // total number of equations: ny + nλ

	// stage: solution and linear solver
	Sol      *Solution   // solution state
	Kb       *la.Triplet // Jacobian == dRdy
	Fb       []float64   // residual == -fb
	Wb       []float64   // workspace
	InitLSol bool        // flag telling that linear solver needs to be initialised prior to any further call
	// contains filtered or unexported fields
}

type DynCoefs struct {
	HHT bool
	// contains filtered or unexported fields
}

type Elem interface {

	// information and initialisation
	Id() int                                           // returns the cell Id
	SetEqs(eqs [][]int, mixedform_eqs []int) (ok bool) // set equations

	// conditions (natural BCs and element's)
	SetEleConds(key string, f fun.Func, extra string) (ok bool) // set element conditions

	// called for each time step
	InterpStarVars(sol *Solution) (ok bool) // interpolate star variables to integration points

	// called for each iteration
	AddToRhs(fb []float64, sol *Solution) (ok bool)                // adds -R to global residual vector fb
	AddToKb(Kb *la.Triplet, sol *Solution, firstIt bool) (ok bool) // adds element K to global Jacobian matrix Kb
	Update(sol *Solution) (ok bool)                                // perform (tangent) update

	// reading and writing of element data
	Encode(enc Encoder) (ok bool) // encodes internal variables
	Decode(dec Decoder) (ok bool) // decodes internal variables

	// output
	OutIpsData() (data []*OutIpData) // returns data from all integration points for output
}

type ElemConnector interface {
	Id() int                                                  // returns the cell Id
	Connect(cid2elem []Elem, c *inp.Cell) (nnzK int, ok bool) // connect multiple elements; e.g.: connect rod/solid elements in Rjoints
}

type ElemIntvars interface {
	Ipoints() (coords [][]float64)                               // returns the real coordinates of integration points [nip][ndim]
	SetIniIvs(sol *Solution, ivs map[string][]float64) (ok bool) // sets initial ivs for given values in sol and ivs map
	BackupIvs(aux bool) (ok bool)                                // create copy of internal variables
	RestoreIvs(aux bool) (ok bool)                               // restore internal variables from copies
	Ureset(sol *Solution) (ok bool)                              // fixes internal variables after u (displacements) have been zeroed
}

type ElemP struct {

	// basic data
	Cid int         // cell/element id
	X   [][]float64 // matrix of nodal coordinates [ndim][nnode]
	Shp *shp.Shape  // shape structure
	Np  int         // total number of unknowns == number of vertices

	// integration points
	IpsElem []*shp.Ipoint // integration points of element
	IpsFace []*shp.Ipoint // integration points corresponding to faces

	// material model
	Mdl *mporous.Model // model

	// problem variables
	Pmap []int // assembly map (location array/element equations)

	// internal variables
	States    []*mporous.State
	StatesBkp []*mporous.State
	StatesAux []*mporous.State

	// gravity
	Gfcn fun.Func // gravity function

	// natural boundary conditions
	NatBcs []*NaturalBc // natural boundary conditions

	Emat     [][]float64 // [nverts][nips] extrapolator matrix
	DoExtrap bool        // do extrapolation of ρl and Cpl => for use with flux and seepage conditions

	// seepage face
	Nf         int   // number of fl variables
	HasSeep    bool  // indicates if this element has seepage faces
	Vid2seepId []int // [nverts] maps local vertex id to index in Fmap
	SeepId2vid []int // [nf] maps seepage face variable id to local vertex id
	Fmap       []int // [nf] map of "fl" variables (seepage face)
	Macaulay   bool  // use discrete ramp function instead of smooth ramp

	Hst   []bool      // [nf] set hydrostatic plmax
	Plmax [][]float64 // [nf][nipsFace] specified plmax (not corrected by multiplier)

	Kpp [][]float64 // [np][np] Kpp := dRpl/dpl consistent tangent matrix
	Kpf [][]float64 // [np][nf] Kpf := dRpl/dfl consistent tangent matrix
	Kfp [][]float64 // [nf][np] Kfp := dRfl/dpl consistent tangent matrix
	Kff [][]float64 // [nf][nf] Kff := dRfl/dfl consistent tangent matrix
	// contains filtered or unexported fields
}

type ElemPhi struct {

	// basic data
	Cid int         // cell/element id
	X   [][]float64 // [ndim][nnode] matrix of nodal coordinates
	Shp *shp.Shape  // shape structure
	Nu  int         // total number of unknowns == number of vertices

	// integration points
	IpsElem []*shp.Ipoint // [nip] integration points of element

	// scratchpad. computed @ each ip
	K [][]float64 // [nu][nu] consistent tangent matrix

	// problem variables
	Umap []int // assembly map (location array/element equations)
	// contains filtered or unexported fields
}

type ElemU struct {

	// basic data
	Cid int         // cell/element id
	X   [][]float64 // matrix of nodal coordinates [ndim][nnode]
	Shp *shp.Shape  // shape structure
	Nu  int         // total number of unknowns

	// variables for dynamics
	Rho  float64  // density of solids
	Cdam float64  // coefficient for damping
	Gfcn fun.Func // gravity function

	// optional data
	UseB      bool    // use B matrix
	Thickness float64 // thickness (for plane-stress)
	Debug     bool    // debugging flag

	// integration points
	IpsElem []*shp.Ipoint // integration points of element
	IpsFace []*shp.Ipoint // integration points corresponding to faces

	// material model and internal variables
	Model    msolid.Model // material model
	MdlSmall msolid.Small // model specialisation for small strains
	MdlLarge msolid.Large // model specialisation for large deformations

	// internal variables
	States    []*msolid.State // [nip] states
	StatesBkp []*msolid.State // [nip] backup states
	StatesAux []*msolid.State // [nip] auxiliary backup states

	// problem variables
	Umap []int // assembly map (location array/element equations)

	// natural boundary conditions
	NatBcs []*NaturalBc

	K [][]float64 // [nu][nu] consistent tangent (stiffness) matrix
	B [][]float64 // [nsig][nu] B matrix for axisymetric case
	D [][]float64 // [nsig][nsig] constitutive consistent tangent matrix

	Δε []float64 // incremental strains leading to updated strains
	// contains filtered or unexported fields
}

type ElemUP struct {

	// auxiliary
	Fconds []*FaceCond // face conditions; e.g. seepage faces
	CtypeU string      // u: cell type
	CtypeP string      // p: cell type

	// underlying elements
	U *ElemU // u-element
	P *ElemP // p-element

	Kup [][]float64 // [nu][np] Kup := dRus/dpl consistent tangent matrix
	Kpu [][]float64 // [np][nu] Kpu := dRpl/dus consistent tangent matrix
	// contains filtered or unexported fields
}

type Encoder interface {
	Encode(e interface{}) error
}

type EssentialBc struct {
	Key   string    // ux, uy, rigid, incsup
	Eqs   []int     // equations
	ValsA []float64 // values for matrix A
	Fcn   fun.Func  // function that implements the "c" in A * y = c
	Inact bool      // inactive
}

type EssentialBcs struct {
	Eq2idx map[int][]int  // maps eq number to indices in BcsTmp
	Bcs    []*EssentialBc // active essential bcs / constraints
	A      la.Triplet     // matrix of coefficients 'A'
	Am     *la.CCMatrix   // compressed form of A matrix

	// temporary
	BcsTmp eqbcpairs // temporary essential bcs / constraints, including inactive ones. maps the first equation number to bcs
}

type FaceCond struct {
	FaceId      int      // msh: cell's face local id
	LocalVerts  []int    // msh: cell's face local vertices ids (sorted)
	GlobalVerts []int    // msh: global vertices ids (sorted)
	Cond        string   // sim: condition; e.g. "qn" or "seepH"
	Func        fun.Func // sim: function to compute boundary condition
	Extra       string   // sim: extra information
}

type GeoLayer struct {
	Tags  []int   // tags of cells within this layer
	Zmin  float64 // coordinate (elevation) at bottom of layer
	Zmax  float64 // coordinate (elevation) at top of layer
	Nodes []*Node // nodes in layer
	Elems []Elem  // elements in layer
	Cl    float64 // liquid compressibility
	RhoS0 float64 // initial density of solids

	K0    float64 // earth-pressure at rest
	Dpl   float64 // liquid pressure added by this layer
	DsigV float64 // absolute value of vertical stress increment added by this layer

	Jac ode.Cb_jac // Jacobian for ode solver
	// contains filtered or unexported fields
}

type GeoLayers []*GeoLayer

type HydroStatic struct {
	Cl float64

	Jac ode.Cb_jac
	// contains filtered or unexported fields
}

type Info struct {

	// essential
	Dofs [][]string        // solution variables PER NODE. ex for 2 nodes: [["ux", "uy", "rz"], ["ux", "uy", "rz"]]
	Y2F  map[string]string // maps "y" keys to "f" keys. ex: "ux" => "fx", "pl" => "ql"

	// internal Dofs; e.g. for mixed formulations
	NintDofs int // number of internal dofs

	// t1 and t2 variables (time-derivatives of first and second order)
	T1vars []string // "pl"
	T2vars []string // "ux", "uy"
}

type NaturalBc struct {
	Key     string   // key such as qn, qn0, ql, seepH, seepP, etc...
	IdxFace int      // local index of face
	Fcn     fun.Func // function callback
	Extra   string   // extra information
}

type Node struct {
	Dofs []*Dof    // degrees-of-freedom == solution variables
	Vert *inp.Vert // pointer to Vertex
}

type OutIpData struct {
	Eid  int                                    // id of element that owns this ip
	X    []float64                              // coordinates
	Calc func(sol *Solution) map[string]float64 // [nkeys] function to calculate secondary values
}

type PtNaturalBc struct {
	Key   string    // key such as fux, fpl, etc...
	Eq    int       // equation
	X     []float64 // location
	Fcn   fun.Func  // function
	Extra string    // extra information
}

type PtNaturalBcs struct {
	Eq2idx map[int]int    // maps eq number to indices in Bcs
	Bcs    []*PtNaturalBc //active boundary conditions such as prescribed forces
}

type RichardsonExtrap struct {

	// variables after big step
	Y_big []float64 // primary variables

	// time loop
	Δt    float64 // time step
	Δtcpy float64 // copy of Δt for divergence control
	// contains filtered or unexported fields
}

type Rjoint struct {

	// basic data
	Edat *inp.ElemData // element data; stored in allocator to be used in Connect
	Cid  int           // cell/element id
	Ny   int           // total number of dofs == rod.Nu + sld.Nu

	// essential
	Rod *Rod            // rod element
	Sld *ElemU          // solid element
	Mdl msolid.RjointM1 // material model

	// shape functions evaluations and extrapolator matrices
	Nmat [][]float64 // [sldNn][rodNn] shape functions of solids @ [N]odes of rod element
	Pmat [][]float64 // [sldNn][rodNp] shape functions of solids @ integration [P]oints of rod element (for Coulomb model)
	Emat [][]float64 // [sldNn][sldNp] solid's extrapolation matrix (for Coulomb model)

	// variables for Coulomb model
	Coulomb bool // use Coulomb model

	// auxiliary variables
	ΔuC [][]float64 // [rodNn][ndim] relative displ. increment of solid @ nodes of rod; Eq (30)
	Δw  []float64   // [ndim] relative velocity; Eq (32)

	// temporary Jacobian matrices. see Eq. (57)
	Krr [][]float64 // [rodNu][rodNu] Eq. (58)
	Krs [][]float64 // [rodNu][sldNu] Eq. (59)
	Ksr [][]float64 // [sldNu][rodNu] Eq. (60)
	Kss [][]float64 // [sldNu][sldNu] Eq. (61)

	// internal values
	States    []*msolid.OnedState // [nip] internal states
	StatesBkp []*msolid.OnedState // [nip] backup internal states
	StatesAux []*msolid.OnedState // [nip] backup internal states
	// contains filtered or unexported fields
}

type Rod struct {

	// basic data
	Cid int         // cell/element id
	X   [][]float64 // matrix of nodal coordinates [ndim][nnode]
	Shp *shp.Shape  // shape structure
	Nu  int         // total number of unknowns == 2 * nsn

	// parameters
	A float64 // cross-sectional area

	// variables for dynamics
	Rho  float64  // density of solids
	Gfcn fun.Func // gravity function

	// integration points
	IpsElem []*shp.Ipoint // integration points of element

	// vectors and matrices
	K   [][]float64 // global K matrix
	M   [][]float64 // global M matrices
	Rus []float64   // residual: Rus = fi - fx

	// problem variables
	Umap []int // assembly map (location array/element equations)

	// material model and internal variables
	Model     msolid.OnedSolid
	States    []*msolid.OnedState
	StatesBkp []*msolid.OnedState
	StatesAux []*msolid.OnedState
	// contains filtered or unexported fields
}

type Solution struct {

	// state
	T      float64   // current time
	Y      []float64 // DOFs (solution variables); e.g. y = {u, p}
	Dydt   []float64 // dy/dt
	D2ydt2 []float64 // d²y/dt²

	// auxiliary
	ΔY  []float64 // total increment (for nonlinear solver)
	Psi []float64 // t1 star vars; e.g. ψ* = β1.p + β2.dpdt
	Zet []float64 // t2 star vars; e.g. ζ* = α1.u + α2.v + α3.a
	Chi []float64 // t2 star vars; e.g. χ* = α4.u + α5.v + α6.a
	L   []float64 // Lagrange multipliers
}

type Summary struct {
	Nproc    int          // number of processors used in last last run; equal to 1 if not distributed
	OutTimes []float64    // [nOutTimes] output times
	Resids   utl.DblSlist // residuals (if Stat is on; includes all stages)
	Dirout   string       // directory where results are stored
	Fnkey    string       // filename key of simulation
}

type T_iteration struct {
	It     int     // iteration number
	ResRel float64 // relative residual
	Resid  float64 // absolute residual
}

type T_results struct {
	Status     string        // status message
	LoadFactor float64       // load factor
	Iterations []T_iteration // iterations data
	Kmats      [][][]float64 // [nele][nu][nu] all stiffness matrices
	Disp       [][]float64   // [nnod][ndim] displacements at nodes
	DispMult   float64       // displacements multiplier
	Note       string        // note about number of integration points
	Sigmas     [][][]float64 // [nele][nip][nsig] all stresses @ all ips 2D:{sx, sy, sxy, sz}
}

type T_results_set []*T_results