Functions

phun(str::String; system_of_units = :SI, base_time_units = :SEC)

Evaluate an expression in physical units.

Inputs: –system_of_units

if system_of_units  ==  :US
   basic assumed units are American Engineering:
   LENGTH = FT, TIME = SEC, MASS = SLUG TEMPERATURE = RAN FORCE = LB
elseif system_of_units  ==  :CGS
   basic assumed units are Centimeter,Gram,Second:
   LENGTH = CM, TIME = SEC, MASS = GM TEMPERATURE = K FORCE = DYNE
elseif system_of_units  ==  :IMPERIAL
   basic assumed units are Imperial:
   LENGTH = FT, TIME = SEC, MASS = SLUG TEMPERATURE = RAN FORCE = LB
otherwise,
   basic assumed units are :SIM (equivalent to :SI, default):
   LENGTH = M , TIME = SEC, MASS = KG   TEMPERATURE = K   FORCE = N

–base_time_units defaults to :SEC

Example

pu = ustring -> phun(ustring; system_of_units = :SIMM)
E1s = 130.0*pu("GPa")

yields 1.3e+5 (in mega Pascal) whereas

130.0*phun("GPa"; system_of_units = :SI)

yields 1.3e+11 (in Pascal)

boundingbox(x::AbstractArray)

Compute the bounding box of the points in x.

x = holds points, one per row.

Returns box = bounding box for 1-D box=[minx,maxx], or for 2-D box=[minx,maxx,miny,maxy], or for 3-D box=[minx,maxx,miny,maxy,minz,maxz]

boxesoverlap(box1::AbstractVector, box2::AbstractVector)

Do the given boxes overlap?

inbox(box::AbstractVector, x::AbstractVector)

Is the given location inside the box?

• box = vector entries arranged as minx,maxx,miny,maxy,minz,maxz.

Note: point on the boundary of the box is counted as being inside.

inflatebox!(box::AbstractVector, inflatevalue::T) where {T}

Inflate the box by the value supplied.

initbox!(box::AbstractVector, x::AbstractVector)

Initialize a bounding box with a single point.

intersectboxes(box1::AbstractVector, box2::AbstractVector)

Compute the intersection of two boxes.

The function returns an empty box (length(b) == 0) if the intersection is empty; otherwise a box is returned.

updatebox!(box::AbstractVector, x::AbstractArray)

Update a box with another location, or create a new box.

If the box does not have the correct dimensions, it is correctly sized.

box = bounding box for 1-D box=[minx,maxx], or for 2-D box=[minx,maxx,miny,maxy], or for 3-D box=[minx,maxx,miny,maxy,minz,maxz] The box is expanded to include the supplied location x. The variable x can hold multiple points in rows.

csmat(self::CSys)

Return coordinate system rotation matrix.

No allocation is involved.

gen_iso_csmat!(csmatout::Matrix{T}, XYZ::Matrix{T}, tangents::Matrix{T}, feid::IT, qpid::IT) where {T, IT}

Compute the coordinate system for an isotropic material using information available by looking at the coordinate curves of isoparametric finite elements.

• XYZ= location in physical coordinates,
• tangents= tangent vector matrix, tangents to the parametric coordinate curves in the element,
• feid= finite element identifier;
• qpid = quadrature point identifier.

The basic assumption here is that the material is isotropic, and therefore the choice of the material directions does not really matter as long as they correspond to the dimensionality of the element. For instance a one-dimensional element (L2 as an example) may be embedded in a three-dimensional space.

This function assumes that it is being called for an mdim-dimensional manifold element, which is embedded in a sdim-dimensional Euclidean space. If mdim == sdim, the coordinate system matrix is the identity; otherwise the local coordinate directions are aligned with the linear subspace defined by the tangent vectors.

updatecsmat!(self::CSys,
    XYZ::Matrix{T},
    tangents::Matrix{T},
    feid::IT1,
    qpid::IT2) where {T, IT1, IT2}

Update the coordinate system orientation matrix.

The coordinate system matrix is updated based upon the location XYZ of the evaluation point, and possibly on the Jacobian matrix tangents within the element in which the coordinate system matrix is evaluated, or perhaps on the identifier feid of the finite element and/or the quadrature point identifier.

After this function returns, the coordinate system matrix can be read in the buffer as self.csmat.

add_b1tdb2!(
    Ke::Matrix{T},
    B1::Matrix{T},
    B2::Matrix{T},
    Jac_w::T,
    D::Matrix{T},
    DB2::Matrix{T},
) where {T}

Add the product (B1'*(D*(Jac_w))*B2), to the matrix Ke.

The matrix Ke is assumed to be suitably initialized: the results of this computation are added. The matrix Ke may be rectangular.

The matrix D may be rectangular.

The matrix Ke is modified. The matrices B1, B2, and D are not modified inside this function. The scratch buffer DB is overwritten during each call of this function.

add_btdb_ut_only!(Ke::Matrix{T}, B::Matrix{T}, Jac_w::T, D::Matrix{T}, DB::Matrix{T}) where {T}

Add the product (B'*(D*(Jac*w[j]))*B), to the matrix Ke.

Only upper triangle is computed; the lower triangle is not touched. (Use complete_lt! to complete the lower triangle, if needed.)

The matrix Ke is assumed to be suitably initialized.

The matrix Ke is modified. The matrices B and D are not modified inside this function. The scratch buffer DB is overwritten during each call of this function.

add_btsigma!(Fe::Vector{T}, B::Matrix{T}, coefficient::T, sigma::Vector{T}) where {T}

Add the product B'*(sigma*coefficient), to the elementwise vector Fe.

The vector Fe is assumed to be suitably initialized.

The vector Fe is modified. The vector sigma is not modified inside this function.

add_gkgt_ut_only!(
    Ke::Matrix{T},
    gradN::Matrix{T},
    Jac_w::T,
    kappa_bar::Matrix{T},
    kappa_bargradNT::Matrix{T},
) where {T}

Add the product gradN*kappa_bar*gradNT*(Jac*w[j]) to the matrix Ke.

Only upper triangle is computed; the lower triangle is not touched. (Use complete_lt! to complete the lower triangle, if needed.)

The matrix Ke is assumed to be suitably initialized.

Upon return, the matrix Ke is updated. The scratch buffer kappa_bargradNT is overwritten during each call of this function. The matrices gradN and kappa_bar are not modified inside this function.

add_mggt_ut_only!(Ke::Matrix{T}, gradN::Matrix{T}, mult) where {T}

Add the product gradN*mult*gradNT to the matrix Ke.

The argument mult is a scalar. Only upper triangle is computed; the lower triangle is not touched. (Use complete_lt! to complete the lower triangle, if needed.)

The matrix Ke is assumed to be suitably initialized.

The matrix Ke is modified. The matrix gradN is not modified inside this function.

add_n1n2t!(Ke::Matrix{T}, N1::Matrix{T}, N2::Matrix{T}, Jac_w_coeff::T) where {T<:Number}

Add the product N1*(N2'*(coeff*(Jac*w(j))), to the matrix Ke.

The matrix Ke is assumed to be suitably initialized. The matrices N1 and N2 have a single column each.

The matrix Ke is modified. The matrix N1 and N2 are not modified inside this function.

add_nnt_ut_only!(Ke::Matrix{T}, N::Matrix{T}, Jac_w_coeff::T) where {T<:Number}

Add the product Nn*(Nn'*(coeff*(Jac*w(j))), to the matrix Ke.

Only the upper triangle is computed; the lower triangle is not touched.

The matrix Ke is assumed to be suitably initialized. The matrix Nn has a single column.

The matrix Ke is modified. The matrix Nn is not modified inside this function.

adjugate3!(B, A)

Compute the adjugate matrix of 3x3 matrix A.

complete_lt!(Ke::Matrix{T}) where {T}

Complete the lower triangle of the elementwise matrix Ke.

The matrix Ke is modified inside this function. The upper-triangle entries are copied across the diagonal to the lower triangle.

detC(::Val{3}, C::Matrix{T})

Compute determinant of 3X3 C.

export_sparse(filnam, M)

Export sparse matrix to a text file.

import_sparse(filnam)

Import sparse matrix from a text file.

jac!(J::Matrix{T}, ecoords::Matrix{T}, gradNparams::Matrix{T}) where {T}

Compute the Jacobian matrix at the quadrature point.

Arguments: J = Jacobian matrix, overwritten inside the function ecoords = matrix of the node coordinates for the element. gradNparams = matrix of basis function gradients

loc!(loc::Matrix{T}, ecoords::Matrix{T}, N::Matrix{T}) where {T}

Compute the location of the quadrature point.

Arguments: loc = matrix of coordinates, overwritten inside the function ecoords = matrix of the node coordinates for the element. N = matrix of basis function values

locjac!(
    loc::Matrix{T},
    J::Matrix{T},
    ecoords::Matrix{T},
    N::Matrix{T},
    gradNparams::Matrix{T},
) where {T}

Compute location and Jacobian matrix at the quadrature point.

Arguments: loc = matrix of coordinates, overwritten inside the function J = Jacobian matrix, overwritten inside the function ecoords = matrix of the node coordinates for the element. N = matrix of basis function values gradNparams = matrix of basis function gradients

matrix_blocked_dd(A, row_nfreedofs, col_nfreedofs = row_nfreedofs)

Extract the "data-data" partition of a matrix.

The matrix is assumed to be composed of four blocks

A = [A_ff A_fd
     A_df A_dd]

Here f stands for free, and d stands for data (i.e. fixed, prescribed, ...). The size of the ff block is row_nfreedofs, col_nfreedofs. Neither one of the blocks is square, unless row_nfreedofs == col_nfreedofs.

When row_nfreedofs == col_nfreedofs, only the number of rows needs to be given.

matrix_blocked_df(A, row_nfreedofs, col_nfreedofs = row_nfreedofs)

Extract the "data-free" partition of a matrix.

The matrix is assumed to be composed of four blocks

A = [A_ff A_fd
     A_df A_dd]

Here f stands for free, and d stands for data (i.e. fixed, prescribed, ...). The size of the ff block is row_nfreedofs, col_nfreedofs. Neither one of the blocks is square, unless row_nfreedofs == col_nfreedofs.

When row_nfreedofs == col_nfreedofs, only the number of rows needs to be given.

matrix_blocked_fd(A, row_nfreedofs, col_nfreedofs = row_nfreedofs)

Extract the "free-data" partition of a matrix.

The matrix is assumed to be composed of four blocks

A = [A_ff A_fd
     A_df A_dd]

Here f stands for free, and d stands for data (i.e. fixed, prescribed, ...). The size of the ff block is row_nfreedofs, col_nfreedofs. Neither one of the blocks is square, unless row_nfreedofs == col_nfreedofs.

When row_nfreedofs == col_nfreedofs, only the number of rows needs to be given.

matrix_blocked_ff(A, row_nfreedofs, col_nfreedofs = row_nfreedofs)

Extract the "free-free" partition of a matrix.

The matrix is assumed to be composed of four blocks

A = [A_ff A_fd
     A_df A_dd]

Here f stands for free, and d stands for data (i.e. fixed, prescribed, ...). The size of the ff block is row_nfreedofs, col_nfreedofs. Neither one of the blocks is square, unless row_nfreedofs == col_nfreedofs.

When row_nfreedofs == col_nfreedofs, only the number of rows needs to be given.

mulCAB!(C, A, B)

Compute the matrix C = A * B

The use of BLAS is purposefully avoided in order to eliminate contentions of multi-threaded execution of the library code with the user-level threads.

Note: See the thread https://discourse.julialang.org/t/ann-loopvectorization/32843/36

mulCAB!(::Val{3}, C, A, B)

Compute the product of 3X3 matrices C = A * B

mulCAB!(C::Vector{T}, A, B::Vector{T})  where {T}

Compute the product C = A * B, where C and B are "vectors".

The use of BLAS is purposefully avoided in order to eliminate contentions of multi-threaded execution of the library code with the user-level threads.

mulCABt!(C, A, B)

Compute the matrix C = A * B'

The use of BLAS is purposefully avoided in order to eliminate contentions of multi-threaded execution of the library code with the user-level threads.

mulCABt!(::Val{3}, C, A, B)

Compute the product of 3X3 matrices C = A * Transpose(B)

mulCAtB!(C, A, B)

Compute the matrix C = A' * B

The use of BLAS is purposefully avoided in order to eliminate contentions of multi-threaded execution of the library code with the user-level threads.

mulCAtB!(::Val{3}, C, A, B)

Compute the product of 3X3 matrices C = Transpose(A) * B

setvectorentries!(a, v = zero(eltype(a)))

Set entries of a long vector to a given constant.

symmetrize!(a)

Make the matrix on input symmetric.

The operation is in-place.

vector_blocked_d(V, nfreedofs)

Extract the "data" part of a vector.

The vector is composed of two blocks

V = [V_f
     V_d]

Here f stands for free, and d stands for data (i.e. fixed, prescribed, ...).

vector_blocked_f(V, nfreedofs)

Extract the "free" part of a vector.

The vector is composed of two blocks

V = [V_f
     V_d]

Here f stands for free, and d stands for data (i.e. fixed, prescribed, ...).

zeros_via_calloc(::Type{T}, dims::Integer...) where {T}

Allocate large array of numbers using calloc.

size(self::DataCache)

Size of the data cache value.

updatenormal!(self::SurfaceNormal, XYZ::Matrix{T}, tangents::Matrix{T}, feid::IT, qpid::IT) where {T, IT}

Update the surface normal vector.

Returns the normal vector (stored in the cache).

updateforce!(self::ForceIntensity, XYZ::Matrix{T}, tangents::Matrix{T}, feid::IT, qpid::IT) where {T<:Number, IT<:Integer}

Update the force intensity vector.

Returns a vector (stored in the cache self.cache).

cross2(theta::AbstractVector{T1}, v::AbstractVector{T2}) where {T1, T2}

Compute the cross product of two vectors in two-space.

cross3!(
    result::AbstractVector{T1},
    theta::AbstractVector{T2},
    v::AbstractVector{T3},
) where {T1, T2, T3}

Compute the cross product of two vectors in three-space in place.

cross3!(
    result::AbstractVector{T1},
    theta::Union{AbstractVector{T2}, Tuple{T2, T2, T2}},
    v::Union{AbstractVector{T3}, Tuple{T3, T3, T3}}
) where {T1, T2, T3}

Compute the cross product of two vectors in three-space in place.

rotmat3!(Rmout::Matrix{T}, a::VT) where {T, VT}

Compute a 3D rotation matrix in-place.

a = array, vector, or tuple with three floating-point numbers

rotmat3(a::VT) where {VT}

Prepare a rotation matrix from a rotation vector

skewmat!(S::Matrix{T}, theta::VT) where {T, VT}

Compute skew-symmetric matrix.

cat(self::T,  other::T) where {T<:AbstractFESet}

Concatenate the connectivities of two FE sets.

count(self::T)::FInt where {T<:AbstractFESet}

Get the number of individual connectivities in the FE set.

eachindex(fes::AbstractFESet)

Create an iterator for elements.

Jacobian(self::ET, J::Matrix{FT}) where {ET<:AbstractFESet0Manifold, FT}

Evaluate the point Jacobian.

• J = Jacobian matrix, columns are tangent to parametric coordinates curves.

Jacobian(self::ET, J::Matrix{FT}) where {ET<:AbstractFESet1Manifold, FT}

Evaluate the curve Jacobian.

• J = Jacobian matrix, columns are tangent to parametric coordinates curves.

Jacobian(self::ET, J::Matrix{FT}) where {ET<:AbstractFESet2Manifold, FT}

Evaluate the curve Jacobian.

• J = Jacobian matrix, columns are tangent to parametric coordinates curves.

Jacobian(self::ET, J::Matrix{FT}) where {ET<:AbstractFESet3Manifold, FT}

Evaluate the volume Jacobian.

J = Jacobian matrix, columns are tangent to parametric coordinates curves.

accepttodelegate(self::T, delegateof) where {T<:AbstractFESet}

Accept to delegate for an object.

bfun(self::ET, param_coords::Vector{T}) where {ET<:AbstractFESet, T}

Compute the values of the basis functions.

Compute the values of the basis functions at a given parametric coordinate. One basis function per row.

bfundpar(self::ET, param_coords::Vector{T}) where {ET<:AbstractFESet, T}

Compute the values of the basis function gradients.

Compute the values of the basis function gradients with respect to the parametric coordinates at a given parametric coordinate. One basis function gradients per row.

boundaryconn(self::T) where {T<:AbstractFESet}

Get boundary connectivity.

boundaryfe(self::T) where {T<:AbstractFESet}

Return the constructor of the type of the boundary finite element.

centroidparametric(self::T) where {T<:AbstractFESet}

Return the parametric coordinates of the centroid of the element.

connasarray(self::AbstractFESet{NODESPERELEM}) where {NODESPERELEM}

Return the connectivity as an array.

Return the connectivity as an integer array (matrix), where the number of rows matches the number of connectivities in the set.

delegateof(self::T) where {T<:AbstractFESet}

Return the object of which the elements set is a delegate.

fromarray!(self::AbstractFESet{NODESPERELEM}, conn::FIntMat) where {NODESPERELEM}

Set the connectivity from an integer array.

gradN!(
    self::AbstractFESet1Manifold,
    gradN::Matrix{FT},
    gradNparams::Matrix{FT},
    redJ::Matrix{FT},
) where {FT}

Compute the gradient of the basis functions with the respect to the "reduced" spatial coordinates.

• gradN= output, matrix of gradients, one per row
• gradNparams= matrix of gradients with respect to parametric coordinates, one per row
• redJ= reduced Jacobian matrix redJ=transpose(Rm)*J

gradN!(
    self::AbstractFESet2Manifold,
    gradN::Matrix{FT},
    gradNparams::Matrix{FT},
    redJ::Matrix{FT},
) where {FT}

Compute the gradient of the basis functions with the respect to the "reduced" spatial coordinates.

• gradN= output, matrix of gradients, one per row
• gradNparams= matrix of gradients with respect to parametric coordinates, one per row
• redJ= reduced Jacobian matrix redJ=transpose(Rm)*J

gradN!(
    self::AbstractFESet3Manifold,
    gradN::Matrix{FT},
    gradNparams::Matrix{FT},
    redJ::Matrix{FT},
) where {FT}

Compute the gradient of the basis functions with the respect to the "reduced" spatial coordinates.

• gradN= output, matrix of gradients, one per row
• gradNparams= matrix of gradients with respect to parametric coordinates, one per row
• redJ= reduced Jacobian matrix redJ=transpose(Rm)*J

inparametric(self::AbstractFESet, param_coords)

Are given parametric coordinates inside the element parametric domain?

Return a Boolean: is the point inside, true or false?

manifdim(me)

Get the manifold dimension.

map2parametric(
    self::ET,
    x::Matrix{FT},
    pt::Vector{FT};
    tolerance = 0.001,
    maxiter = 5,
) where {ET<:AbstractFESet, FT}

Map a spatial location to parametric coordinates.

• x=array of spatial coordinates of the nodes, size(x) = nbfuns x dim,
• c= spatial location
• tolerance = tolerance in parametric coordinates; default is 0.001.

Return

• success = Boolean flag, true if successful, false otherwise.
• pc = Returns a row array of parametric coordinates if the solution was successful, otherwise NaN are returned.

nodesperelem(fes::AbstractFESet{NODESPERELEM}) where {NODESPERELEM}

Provide the number of nodes per element.

nodesperelem(::Type{T}) where {NODESPERELEM, T<:AbstractFESet{NODESPERELEM}}

Provide the number of nodes per element for a given type.

setlabel!(self::ET, val::IT) where {ET<:AbstractFESet, IT}

Set the label of the entire finite elements set.

All elements are labeled with this number.

setlabel!(self::ET, val::Vector{IT}) where {ET<:AbstractFESet, IT}

Set the labels of individual elements.

subset(self::T, L) where {T<:AbstractFESet}

Extract a subset of the finite elements from the given finite element set.

• L: an integer vector, tuple, or a range.

updateconn!(self::ET, newids::Vector{IT}) where {ET<:AbstractFESet, IT}

Update the connectivity after the IDs of nodes changed.

newids= new node IDs. Note that indexes in the conn array "point" into the newids array. After the connectivity was updated this will no longer be true!

count(self::FENodeSet)

Get the number of finite element nodes in the node set.

eachindex(fens::FENodeSet)

Create the finite element node iterator.

spacedim(self::FENodeSet)

Number of dimensions of the space in which the node lives, 1, 2, or 3.

xyz3(self::FENodeSet)

Get the 3-D coordinate that define the location of the node. Even if the nodes were specified in lower dimension (1-D, 2-D) this function returns a 3-D coordinate by padding with zeros.

connectedelems(fes::AbstractFESet, node_list::FIntVec)

Extract the list of numbers for the fes that are connected to given nodes.

connectednodes(fes::AbstractFESet)

Extract the node numbers of the nodes connected by given finite elements.

Extract the list of unique node numbers for the nodes that are connected by the finite element set fes. Note that it is assumed that all the FEs are of the same type (the same number of connected nodes by each cell).

findunconnnodes(fens::FENodeSet, fes::AbstractFESet)

Find nodes that are not connected to any finite element.

connected = array is returned which is for the node k either true (node k is connected), or false (node k is not connected).

Let us say there are nodes not connected to any finite element that you would like to remove from the mesh: here is how that would be accomplished.

selectelem(fens::FENodeSet, fes::T; kwargs...) where {T<:AbstractFESet}

Select finite elements.

Arguments

• fens = finite element node set
• fes = finite element set
• kwargs = keyword arguments to specify the selection criteria

Selection criteria

facing

Select all "boundary" elements that "face" a certain direction.

exteriorbfl = selectelem(fens, bdryfes, facing=true, direction=[1.0, 1.0, 0.0]);

or

exteriorbfl = selectelem(fens, bdryfes, facing=true, direction=dout, dotmin = 0.99);

where

function dout(xyz)
    return xyz/norm(xyz)
end

and xyz is the location of the centroid of a boundary element. Here the finite element is considered "facing" in the given direction if the dot product of its normal and the direction vector is greater than dotmin. The default value for dotmin is 0.01 (this corresponds to almost 90 degrees between the normal to the finite element and the given direction).

This selection method makes sense only for elements that are surface-like (i. e. for boundary mmeshes).

label

Select elements based on their label.

rl1 = selectelem(fens, fes, label=1)

box, distance

Select elements based on some criteria that their nodes satisfy. See the function selectnode().

Example: Select all elements whose nodes are closer than R+inflate from the point from.

linner = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0],
  inflate = tolerance)

Example:

exteriorbfl = selectelem(fens, bdryfes,
   box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);

withnodes

Select elements whose nodes are in a given list of node numbers.

Example:

l = selectelem(fens, fes, withnodes = [13, 14])

flood

Select all FEs connected together, starting from a given node. Connections through a vertex (node) are sufficient.

Example: Select all FEs connected together (Starting from node 13):

l = selectelem(fens, fes, flood = true, startnode = 13)

Optional keyword arguments

Should we consider the element only if all its nodes are in?

• allin = Boolean: if true, then all nodes of an element must satisfy

the criterion; otherwise one is enough.

Output

felist = list of finite elements from the set that satisfy the criteria

selectnode(fens::FENodeSet; kwargs...)

Select nodes using some criterion.

Arguments

• v = array of locations, one location per row
• kwargs = pairs of keyword argument/value

Selection criteria

box

nLx = vselect(fens.xyz, box = [0.0 Lx  0.0 0.0 0.0 0.0], inflate = Lx/1.0e5)

The keyword 'inflate' may be used to increase or decrease the extent of the box (or the distance) to make sure some nodes which would be on the boundary are either excluded or included.

distance

list = selectnode(fens.xyz; distance=1.0+0.1/2^nref, from=[0. 0.],
        inflate=tolerance);

Find all nodes within a certain distance from a given point.

plane

candidates = selectnode(fens; plane = [0.0 0.0 1.0 0.0], thickness = h/1000)

The keyword plane defines the plane by its normal (the first two or three numbers) and its distance from the origin (the last number). Nodes are selected they lie on the plane, or near the plane within the distance thickness from the plane. The normal is assumed to be of unit length, if it isn't apply as such, it will be normalized internally.

nearestto

nh = selectnode(fens; nearestto = [R+Ro/2, 0.0, 0.0] )

Find the node nearest to the location given.

farthestfrom

nh = selectnode(fens; farthestfrom = [R+Ro/2, 0.0, 0.0] )

Find the node farthest from the location given.

vselect(v::Matrix{T}; kwargs...) where {T<:Number}

Select locations (vertices) from the array based on some criterion.

See the function selectnode() for examples of the criteria that can be used to search vertices.

copyto!(DEST::F,  SRC::F) where {F<:AbstractField}

Copy data from one field to another.

anyfixedvaluenz(self::F, conn::CC) where {F<:AbstractField, CC}

Is any degree of freedom fixed (prescribed) to be non-zero?

applyebc!(self::F) where {F<:AbstractField}

Apply EBCs (essential boundary conditions).

fixeddofs(self::F) where {F<:AbstractField}

Return range corresponding to the fixed degrees of freedom.

freedofs(self::F) where {F<:AbstractField}

Return range corresponding to the free degrees of freedom.

gatherdofnums!(self::F, dest::A, conn::CC) where {F<:AbstractField, A, CC}

Gather dofnums from the field.

The order is: for each node in the connectivity, copy into the buffer all the degrees of freedom for that node, then the next node and so on.

gatherfixedvalues_asmat!(
    self::F,
    dest::AbstractArray{T,2},
    conn::CC,
) where {F<:AbstractField, T, CC}

Gather FIXED values from the field into a two-dimensional array.

The order is: for each node in the connectivity, copy into the corresponding row of the buffer all the degrees of freedom, then the next node into the next row and so on. If a degree of freedom is NOT fixed, the corresponding entry is set to zero.

dest = destination buffer: overwritten inside, must be preallocated in the correct size

gatherfixedvalues_asvec!(
    self::F,
    dest::AbstractArray{T,1},
    conn::CC,
) where {F<:AbstractField, T, CC}

Gather FIXED values from the field into a vector.

The order is: for each node in the connectivity, copy into the buffer all the fixed degrees of freedom, then the next node and so on. If a degree of freedom is NOT fixed, the corresponding entry is set to zero.

dest = destination buffer: overwritten inside, must be preallocated in the correct size

gathersysvec!(self::F,
    vec::Vector{T}, which = :f) where {F<:AbstractField, T}

Gather values from the field for the system vector.

Arguments

• self: field;
• which:
  • :f - Collect a vector that includes the free degrees of freedom (default);
  • :d - Collect a vector that includes the fixed (data) degrees of freedom;
  • :a - Collect a vector that includes all the degrees of freedom, free and fixed.

The system vector consists of two parts: the first part, from 1 to nfreedofs (self) are the free degrees of freedom, the second part from nfreedofs (self)+1 to nalldofs(self) are the fixed degrees of freedom.

This function gathers either the entire vector, or one of the parts. The length of the supplied buffer vec must be correct, either nfreedofs (self), nalldofs(self), or nfixeddofs(self).

gathersysvec(self::F, which = :f) where {F<:AbstractField}

Gather values from the field for the system vector.

Arguments

• self: field;
• which:
  • :f - Collect a vector that includes the free degrees of freedom (default);
  • :d - Collect a vector that includes the fixed (data) degrees of freedom;
  • :a - Collect a vector that includes all the degrees of freedom, free and fixed.

The system vector consists of two parts: the first part, from 1 to nfreedofs (self) are the free degrees of freedom, the second part from nfreedofs (self)+1 to nalldofs(self) are the fixed degrees of freedom.

This function returns either the entire vector, or one of the parts.

gathervalues_asmat!(
    self::F,
    dest::AbstractArray{T,2},
    conn::CC,
) where {F<:AbstractField, T, CC}

Gather values from the field into a two-dimensional array.

The order is: for each node in the connectivity, copy into the corresponding row of the buffer all the degrees of freedom, then the next node into the next row and so on.

dest = destination buffer: overwritten inside, must be preallocated in the correct size

gathervalues_asvec!(
    self::F,
    dest::AbstractArray{T,1},
    conn::CC,
) where {F<:AbstractField, T, CC}

Gather values from the field into a vector.

The order is: for each node in the connectivity, copy into the buffer all the degrees of freedom, then the next node and so on.

dest = destination buffer: overwritten inside, must be preallocated in the correct size

incrscattersysvec!(self::F, vec::AbstractVector{T}) where {F<:AbstractField, T<:Number}

Increment values of the field by scattering a system vector.

The vector may be either for just the free degrees of freedom, or for all the degrees of freedom.

nalldofs(self::F)

Return to number of ALL degrees of freedom (total number of degrees of freedom, which is equal to the number of degrees of freedom per entity times the number of entities).

ndofs(self::F)

Dimension of the degree of freedom parameters (i. e. how many degrees of freedom per entity).

nents(self::F)

Number of nodes associated with the field.

nfixeddofs(self::F)

Return to number of FIXED degrees of freedom (known, data).

nfreedofs(self::F)

Return to number of FREE degrees of freedom.

numberdofs!(self::F) where {F<:AbstractField}

Number the degrees of freedom.

The free components in the field are numbered consecutively, then all the fixed components are numbered, again consecutively.

No effort is made to optimize the numbering in any way. If you'd like to optimize the numbering of the degrees of freedom, use a form that sets the permutation of the degrees of freedom, or the permutation of the nodes.

numberdofs!(self::F, entperm) where {F<:AbstractField}

Number the degrees of freedom.

The free components in the field are numbered consecutively, then all the fixed components are numbered, again consecutively.

The sequence of the entities is given by the entperm permutation (array or range).

prescribeddofs(uebc::F1, u::F2) where {F1<:AbstractField,  F2<:AbstractField}

Find which degrees of freedom are prescribed. uebc = field which defines the constraints (is the dof fixed and to which value), u = field which does not have the constraints applied, and serves as the source of equation numbers; uebc and u may be one and the same field.

scattersysvec!(self::F, vec::AbstractVector{T}) where {F<:AbstractField, T<:Number}

Scatter values to the field from a system vector.

The vector may be either for just the free degrees of freedom, or for all the degrees of freedom.

setebc!(self::F)

Set the EBCs (essential boundary conditions).

All essential boundary conditions are CLEARED.

Note: Any call to setebc!() potentially changes the current assignment which degrees of freedom are free and which are fixed and therefore is presumed to invalidate the current degree-of-freedom numbering. In such a case this method sets _nfreedofs = 0; and dofnums=0.

setebc!(self::F, fenids::AbstractVector{IT})  where {IT<:Integer}

Set the EBCs (essential boundary conditions).

Suppress all degrees of freedom at the given nodes.

fenids - array of N node identifiers

Note: Any call to setebc!() potentially changes the current assignment which degrees of freedom are free and which are fixed and therefore is presumed to invalidate the current degree-of-freedom numbering. In such a case this method sets _nfreedofs = 0; and dofnums=0.

setebc!(self::F, fenid::IT) where {IT<:Integer}

Set the EBCs (essential boundary conditions).

Suppress all degrees of freedom at the given node.

fenid - One integer as a node identifier

Note: Any call to setebc!() potentially changes the current assignment which degrees of freedom are free and which are fixed and therefore is presumed to invalidate the current degree-of-freedom numbering. In such a case this method sets _nfreedofs = 0; and dofnums=0.

setebc!(
    self::F,
    fenids::AbstractVector{IT},
    is_fixed::Bool,
    comp::AbstractVector{IT},
    val::T = 0.0,
) where {T<:Number, IT<:Integer}

Set the EBCs (essential boundary conditions).

fenids = array of N node identifiers comp = integer vector, which degree of freedom (component), val = scalar of type T, default is 0.0

Note: Any call to setebc!() potentially changes the current assignment which degrees of freedom are free and which are fixed and therefore is presumed to invalidate the current degree-of-freedom numbering. In such a case this method sets _nfreedofs = 0; and dofnums=0.

setebc!(
    self::F,
    fenids::AbstractVector{IT},
    is_fixed::Bool,
    comp::IT,
    val::AbstractVector{T},
) where {T<:Number, IT<:Integer}

Set the EBCs (essential boundary conditions).

fenids - array of N node identifiers is_fixed = scaler Boolean: are the degrees of freedom being fixed (true) or released (false), comp = integer, which degree of freedom (component), val = array of N values of type T

Note: Any call to setebc!() potentially changes the current assignment which degrees of freedom are free and which are fixed and therefore is presumed to invalidate the current degree-of-freedom numbering. In such a case this method sets _nfreedofs = 0; and dofnums=0.

setebc!(
    self::F,
    fenids::AbstractVector{IT},
    is_fixed::Bool,
    comp::IT,
    val::T = 0.0,
) where {T<:Number, IT<:Integer}

Set the EBCs (essential boundary conditions).

fenids - array of N node identifiers is_fixed = scaler Boolean: are the degrees of freedom being fixed (true) or released (false), comp = integer, which degree of freedom (component), val = scalar of type T

Note: Any call to setebc!() potentially changes the current assignment which degrees of freedom are free and which are fixed and therefore is presumed to invalidate the current degree-of-freedom numbering. In such a case this method sets _nfreedofs = 0; and dofnums=0.

setebc!(
    self::F,
    fenids::AbstractVector{IT},
    comp::IT,
    val::AbstractVector{T},
) where {T<:Number, IT<:Integer}

Set the EBCs (essential boundary conditions).

fenids = array of N node identifiers comp = integer, which degree of freedom (component), val = array of N values of type T

Note: Any call to setebc!() potentially changes the current assignment which degrees of freedom are free and which are fixed and therefore is presumed to invalidate the current degree-of-freedom numbering. In such a case this method sets _nfreedofs = 0; and dofnums=0.

setebc!(
    self::F,
    fenids::AbstractVector{IT},
    comp::IT,
    val::T = 0.0,
) where {T<:Number, IT<:Integer}

Set the EBCs (essential boundary conditions).

fenids = array of N node identifiers comp = integer, which degree of freedom (component), val = scalar of type T

Note: Any call to setebc!() potentially changes the current assignment which degrees of freedom are free and which are fixed and therefore is presumed to invalidate the current degree-of-freedom numbering. In such a case this method sets _nfreedofs = 0; and dofnums=0.

setebc!(
    self::F,
    fenid::IT,
    is_fixed::Bool,
    comp::IT,
    val::T,
) where {T<:Number, IT<:Integer}

Set the EBCs (essential boundary conditions).

fenids - array of N node identifiers is_fixed = scaler Boolean: are the degrees of freedom being fixed (true) or released (false), comp = integer, which degree of freedom (component), val = array of N values of type T

Note: Any call to setebc!() potentially changes the current assignment which degrees of freedom are free and which are fixed and therefore is presumed to invalidate the current degree-of-freedom numbering. In such a case this method sets _nfreedofs = 0; and dofnums=0.

wipe!(self::F) where {F<:AbstractField}

Wipe all the data from the field.

This includes values, prescribed values, degree of freedom numbers, and "is fixed" flags. The number of free degrees of freedom is set to zero.

nnodes(self::NodalField)

Provide the number of nodes in the nodal field.

nelems(self::ElementalField)

Provide the number of elements in the elemental field.

Jacobiancurve(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
) where {MT<:AbstractFESet0Manifold, CC, T<:Number}

Evaluate the curve Jacobian.

• J = Jacobian matrix
• loc = location of the quadrature point in physical coordinates,
• conn = connectivity of the element,
• N = matrix of basis function values at the quadrature point.

Jacobiancurve(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
) where {MT<:AbstractFESet1Manifold, CC, T<:Number}

Evaluate the curve Jacobian.

• J = Jacobian matrix
• loc = location of the quadrature point in physical coordinates,
• conn = connectivity of the element,
• N = matrix of basis function values at the quadrature point.

Jacobianmdim(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
    m::IT,
) where {MT<:AbstractFESet0Manifold, CC, T<:Number, IT}

Evaluate the manifold Jacobian for an m-dimensional manifold.

For an 0-dimensional finite element, the manifold Jacobian is for

• m=0: +1
• m=1: Jacobiancurve
• m=2: Jacobiansurface
• m=3: Jacobianvolume

Jacobianmdim(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
    m::IT,
) where {MT<:AbstractFESet1Manifold, CC, T<:Number, IT}

Evaluate the manifold Jacobian for an m-dimensional manifold.

For an 1-dimensional finite element, the manifold Jacobian is for

• m=1: Jacobiancurve
• m=2: Jacobiansurface
• m=3: Jacobianvolume

Jacobianmdim(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
    m::IT,
) where {MT<:AbstractFESet2Manifold, CC, T<:Number, IT}

Evaluate the manifold Jacobian for an m-dimensional manifold.

For an 2-dimensional finite element, the manifold Jacobian is for

• m=2: Jacobiansurface
• m=3: Jacobianvolume

Jacobianmdim(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
    m::IT,
) where {MT<:AbstractFESet3Manifold, CC, T<:Number, IT}

Evaluate the manifold Jacobian for an m-dimensional manifold.

For an 3-dimensional cell, the manifold Jacobian is

• m=3: Jacobianvolume

Jacobianpoint(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
) where {MT<:AbstractFESet0Manifold, CC, T<:Number}

Evaluate the point Jacobian.

• J = Jacobian matrix
• loc = location of the quadrature point in physical coordinates,
• conn = connectivity of the element,
• N = matrix of basis function values at the quadrature point.

Jacobiansurface(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
) where {MT<:AbstractFESet0Manifold, CC, T<:Number}

Evaluate the surface Jacobian.

For the zero-dimensional cell, the surface Jacobian is (i) the product of the point Jacobian and the other dimension (units of length squared); or, when used as axially symmetric (ii) the product of the point Jacobian and the circumference of the circle through the point loc times the other dimension (units of length).

• J = Jacobian matrix
• loc = location of the quadrature point in physical coordinates,
• conn = connectivity of the element,
• N = matrix of basis function values at the quadrature point.

Jacobiansurface(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
) where {MT<:AbstractFESet1Manifold, CC, T<:Number}

Evaluate the surface Jacobian.

For the one-dimensional cell, the surface Jacobian is (i) the product of the curve Jacobian and the other dimension (units of length); or, when used as axially symmetric (ii) the product of the curve Jacobian and the circumference of the circle through the point loc.

• J = Jacobian matrix
• loc = location of the quadrature point in physical coordinates,
• conn = connectivity of the element,
• N = matrix of basis function values at the quadrature point.

Jacobiansurface(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
) where {MT<:AbstractFESet2Manifold, CC, T<:Number}

Evaluate the surface Jacobian.

• J = Jacobian matrix
• loc = location of the quadrature point in physical coordinates,
• conn = connectivity of the element,
• N = matrix of basis function values at the quadrature point.

Jacobianvolume(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
) where {MT<:AbstractFESet0Manifold, CC, T<:Number}

Evaluate the volume Jacobian.

For the zero-dimensional cell, the volume Jacobian is (i) the product of the point Jacobian and the other dimension (units of length cubed); or, when used as axially symmetric (ii) the product of the point Jacobian and the circumference of the circle through the point loc and the other dimension (units of length squared).

• J = Jacobian matrix
• loc = location of the quadrature point in physical coordinates,
• conn = connectivity of the element,
• N = matrix of basis function values at the quadrature point.

Jacobianvolume(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
) where {MT<:AbstractFESet1Manifold, CC, T<:Number}

Evaluate the volume Jacobian.

For the one-dimensional cell, the volume Jacobian is (i) the product of the curve Jacobian and the other dimension (units of length squared); or, when used as axially symmetric (ii) the product of the curve Jacobian and the circumference of the circle through the point loc and the other dimension (units of length).

• J = Jacobian matrix
• loc = location of the quadrature point in physical coordinates,
• conn = connectivity of the element,
• N = matrix of basis function values at the quadrature point.

Jacobianvolume(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
) where {MT<:AbstractFESet2Manifold, CC, T<:Number}

Evaluate the volume Jacobian.

For the two-dimensional cell, the volume Jacobian is (i) the product of the surface Jacobian and the other dimension (units of length); or, when used as axially symmetric (ii) the product of the surface Jacobian and the circumference of the circle through the point loc (units of length).

• J = Jacobian matrix
• loc = location of the quadrature point in physical coordinates,
• conn = connectivity of the element,
• N = matrix of basis function values at the quadrature point.

Jacobianvolume(
    self::IntegDomain{MT},
    J::Matrix{T},
    loc::Matrix{T},
    conn::CC,
    N::Matrix{T},
) where {MT<:AbstractFESet3Manifold, CC, T<:Number}

Evaluate the volume Jacobian.

• J = Jacobian matrix
• loc = location of the quadrature point in physical coordinates,
• conn = connectivity of the element,
• N = matrix of basis function values at the quadrature point.

integrationdata(self::IntegDomain)

Calculate the data needed for numerical quadrature for the integration rule stored by the integration domain.

integrationdata(
    self::IntegDomain,
    integration_rule::IR,
) where {IR<:AbstractIntegRule}

Calculate the data needed for a given numerical quadrature rule.

For given integration domain, compute the quantities needed for numerical integration. The integration rule does not necessarily have to be the one associated originally with the integration domain.

Return

npts, Ns, gradNparams, w, pc = number of quadrature points, arrays of basis function values at the quadrature points, arrays of gradients of basis functions with respect to the parametric coordinates, array of weights and array of locations of the quadrature points.

otherdimensionunity(loc::Matrix{T}, conn::CC, N::Matrix{T}) where {CC, T<:Number}

Evaluate the other dimension: default is 1.0.

eltype(a::A) where {A <: AbstractSysmatAssembler}

What is the type of the matrix buffer entries?

assemble!(self::SysmatAssemblerFFBlock,
    mat::MBT,
    dofnums_row::CIT,
    dofnums_col::CIT) where {MBT, CIT}

Assemble a matrix, v

assemble!(
    self::SysmatAssemblerSparse,
    mat::MT,
    dofnums_row::IT,
    dofnums_col::IT,
) where {MT, IT}

Assemble a rectangular matrix.

assemble!(
    self::SysmatAssemblerSparseDiag,
    mat::MT,
    dofnums_row::IV,
    dofnums_col::IV,
) where {MT, IV}

Assemble a square symmetric diagonal matrix.

• dofnums = the row degree of freedom numbers, the column degree of freedom

number input is ignored (the row and column numbers are assumed to be the same).

• mat = diagonal square matrix

assemble!(
    self::SysmatAssemblerSparseHRZLumpingSymm,
    mat::MT,
    dofnums::IV,
    ignore::IV,
) where {MT, IV}

Assemble a HRZ-lumped square symmetric matrix.

Assembly of a HRZ-lumped square symmetric matrix. The method assembles the scaled diagonal of the square symmetric matrix using the two vectors of equation numbers for the rows and columns.

assemble!(
    self::SysmatAssemblerSparseSymm,
    mat::MT,
    dofnums::IT,
    ignore
) where {MT, IT}

Assemble a square symmetric matrix.

dofnums are the row degree of freedom numbers, the column degree of freedom number input is ignored (the row and column numbers are assumed to be the same).

assemble!(self::SysvecAssembler{T}, vec::MV,
  dofnums::D) where {T<:Number, MV<:AbstractArray{T}, D<:AbstractArray{FInt}}

Assemble an elementwise vector.

The method assembles a column element vector using the vector of degree of freedom numbers for the rows.

assemble!(self::SysvecAssembler{T}, vec::MV,
  dofnums::D) where {T<:Number, MV<:AbstractArray{T}, D<:AbstractArray{FInt}}

Assemble an elementwise vector.

The method assembles a column element vector using the vector of degree of freedom numbers for the rows.

assemble!(self::SysvecAssemblerFBlock,
    vec::MV,
    dofnums::IV) where {MV, IV}

Assemble an elementwise vector.

expectedntriples(
    a::A,
    elem_mat_nrows::IT,
    elem_mat_ncols::IT,
    n_elem_mats::IT,
) where {A<:AbstractSysmatAssembler, IT}

How many triples (I, J, V) does the assembler expect to store?

Default is: the product of the size of the element matrices times the number of matrices.

makematrix!(self::SysmatAssemblerFFBlock)

Make an assembled matrix. Delegate the construction of the matrix to the wrapped assembler. Then extract the left upper corner block of the matrix(the free-free matrix).

makematrix!(self::SysmatAssemblerSparseDiag)

Make a sparse symmetric square diagonal matrix.

makematrix!(self::SysmatAssemblerSparseHRZLumpingSymm)

Make a sparse HRZ-lumped symmetric square matrix.

makematrix!(self::SysmatAssemblerSparseSymm)

Make a sparse symmetric square matrix.

makematrix!(self::SysmatAssemblerSparse)

Make a sparse matrix.

The sparse matrix is returned.

makevector!(self::SysvecAssembler)

Make the global vector.

makevector!(self::SysvecAssemblerFBlock)

Make the global "free" vector.

makevector!(self::SysvecAssembler)

Make the global vector.

startassembly!(self::SysvecAssembler{T}, ndofs_row::FInt) where {T<:Number}

Start assembly.

The method makes the buffer for the vector assembly. It must be called before the first call to the method assemble.

• elem_mat_nmatrices = number of element matrices expected to be processed during the assembly.
• ndofs_row= Total number of degrees of freedom.

Returns

• self: the modified assembler.

startassembly!(self::SysmatAssemblerFFBlock,
    elem_mat_nrows::IT,
    elem_mat_ncols::IT,
    n_elem_mats::IT,
    row_nalldofs::IT,
    col_nalldofs::IT;
    force_init = false) where {IT <: Integer}

Start assembly, delegate to the wrapped assembler.

startassembly!(self::SysvecAssembler, ndofs_row)

Start assembly.

The method makes the buffer for the vector assembly. It must be called before the first call to the method assemble.

ndofs_row= Total number of degrees of freedom.

Returns

• self: the modified assembler.

startassembly!(self::SysvecAssemblerFBlock, row_nalldofs::IT) where {IT <: Integer}

Start assembly.

startassembly!(self::SysmatAssemblerSparseDiag{T},
    elem_mat_nrows::IT,
    elem_mat_ncols::IT,
    n_elem_mats::IT,
    row_nalldofs::IT,
    col_nalldofs::IT;
    force_init = false
    ) where {T, IT<:Integer}

Start the assembly of a symmetric square diagonal matrix.

The method makes buffers for matrix assembly. It must be called before the first call to the method assemble!.

Arguments

• elem_mat_nrows = row dimension of the element matrix;
• elem_mat_ncols = column dimension of the element matrix;
• n_elem_mats = number of element matrices;
• row_nalldofs= The total number of rows as a tuple;
• col_nalldofs= The total number of columns as a tuple.

The values stored in the buffers are initially undefined!

Returns

• self: the modified assembler.

startassembly!(self::SysmatAssemblerSparseHRZLumpingSymm{T},
        elem_mat_nrows::IT,
        elem_mat_ncols::IT,
        n_elem_mats::IT,
        row_nalldofs::IT,
        col_nalldofs::IT;
        force_init = false
        ) where {T, IT<:Integer}

Start the assembly of a symmetric lumped diagonal square global matrix.

The method makes buffers for matrix assembly. It must be called before the first call to the method assemble!.

Arguments

• elem_mat_nrows = row dimension of the element matrix;
• elem_mat_ncols = column dimension of the element matrix;
• n_elem_mats = number of element matrices;
• row_nalldofs= The total number of rows as a tuple;
• col_nalldofs= The total number of columns as a tuple.

The values stored in the buffers are initially undefined!

Returns

• self: the modified assembler.

startassembly!(self::SysmatAssemblerSparseSymm{T},
    elem_mat_nrows::IT,
    elem_mat_ncols::IT,
    n_elem_mats::IT,
    row_nalldofs::IT,
    col_nalldofs::IT;
    force_init = false
    ) where {T, IT<:Integer}

Start the assembly of a global matrix.

The method makes buffers for matrix assembly. It must be called before the first call to the method assemble!.

Arguments

• elem_mat_nrows = row dimension of the element matrix;
• elem_mat_ncols = column dimension of the element matrix;
• n_elem_mats = number of element matrices;
• row_nalldofs= The total number of rows as a tuple;
• col_nalldofs= The total number of columns as a tuple.

If the buffer_pointer field of the assembler is 0, which is the case after that assembler was created, the buffers are resized appropriately given the dimensions on input. Otherwise, the buffers are left completely untouched. The buffer_pointer is set to the beginning of the buffer, namely 1.

Returns

• self: the modified assembler.

startassembly!(self::SysmatAssemblerSparse{T},
    elem_mat_nrows::IT,
    elem_mat_ncols::IT,
    n_elem_mats::IT,
    row_nalldofs::IT,
    col_nalldofs::IT;
    force_init = false
    ) where {T, IT<:Integer}

Start the assembly of a global matrix.

The method makes buffers for matrix assembly. It must be called before the first call to the method assemble!.

Arguments

• elem_mat_nrows = row dimension of the element matrix;
• elem_mat_ncols = column dimension of the element matrix;
• n_elem_mats = number of element matrices;
• row_nalldofs= The total number of rows as a tuple;
• col_nalldofs= The total number of columns as a tuple.

If the buffer_pointer field of the assembler is 0, which is the case after that assembler was created, the buffers are resized appropriately given the dimensions on input. Otherwise, the buffers are left completely untouched. The buffer_pointer is set to the beginning of the buffer, namely 1.

Returns

• self: the modified assembler.

L2block(Length::T, nL::IT) where {T<:Number, IT<:Integer}

Mesh of a 1-D block of L2 finite elements.

L2blockx(xs::Vector{T}) where {T<:Number}

Graded mesh of a 1-D block, L2 finite elements.

L3blockx(xs::Vector{T}) where {T<:Number}

Graded mesh of a 1-D block, L2 finite elements.

Q4toT3(fens::FENodeSet, fes::FESetQ4, orientation::Symbol=:default)

Convert a mesh of quadrilateral Q4's to two T3 triangles each.

T3annulus(rin::T, rex::T, nr::IT, nc::IT, angl::T, orientation::Symbol=:a)

Mesh of an annulus segment.

Mesh of an annulus segment, centered at the origin, with internal radius rin, and external radius rex, and development angle angl (in radians). Divided into elements: nr, nc in the radial and circumferential direction respectively.

T3block(Length::T, Width::T, nL::IT, nW::IT, orientation::Symbol = :a
    ) where {T<:Number, IT<:Integer}

T3 mesh of a rectangle.

T3blockx(xs::VecOrMat{T}, ys::VecOrMat{T}, orientation::Symbol = :a
    ) where {T<:Number}

T3 mesh of a rectangle.

T3circlen(radius::T, nperradius::IT) where {T<:Number, IT<:Integer}

Mesh of a quarter circle with a given number of elements per radius.

The parameter nperradius should be an even number; if that isn't so is adjusted to by adding one.

T3circleseg(
    angle::T,
    radius::T,
    ncircumferentially::IT,
    nperradius::IT,
    orientation::Symbol = :a,
) where {T<:Number, IT<:Integer}

Mesh of a segment of a circle.

The subtended angle is angle in radians. The orientation: refer to T3block.

T3refine(fens::FENodeSet,fes::FESetT3)

Refine a mesh of 3-node triangles by quadrisection.

T3toT6(fens::FENodeSet, fes::FESetT3)

Convert a mesh of triangle T3 (three-node) to triangle T6.

T6annulus(
    rin::T,
    rex::T,
    nr::IT,
    nc::IT,
    angl::T,
    orientation::Symbol = :a,
) where {T<:Number, IT<:Integer}

Mesh of an annulus segment.

Mesh of an annulus segment, centered at the origin, with internal radius rin, and external radius rex, and development angle angl (in radians). Divided into elements: nr, nc in the radial and circumferential direction respectively.

T6block(Length::T, Width::T, nL::IT, nW::IT, orientation::Symbol = :a
    ) where {T<:Number, IT<:Integer}

Mesh of a rectangle of T6 elements.

T6blockx(xs::VecOrMat{T}, ys::VecOrMat{T}, orientation::Symbol = :a
    ) where {T<:Number}

Graded mesh of a 2-D block of T6 finite elements.

Q4annulus(rin::T, rex::T, nr::IT, nc::IT, Angl::T) where {T<:Number, IT<:Integer}

Mesh of an annulus segment.

Mesh of an annulus segment, centered at the origin, with internal radius rin, and external radius rex, and development angle Angl (in radians). Divided into elements: nr, nc in the radial and circumferential direction respectively.

Q4block(Length::T, Width::T, nL::IT, nW::IT) where {T<:Number, IT<:Integer}

Mesh of a rectangle, Q4 elements.

Divided into elements: nL, nW in the first, second (x,y).

Q4blockx(xs::Vector{T}, ys::Vector{T})

Graded mesh of a rectangle, Q4 finite elements.

Mesh of a 2-D block, Q4 finite elements. The nodes are located at the Cartesian product of the two intervals on the input. This allows for construction of graded meshes.

xs,ys - Locations of the individual planes of nodes.

Q4circlen(radius::T, nperradius::IT) where {T<:Number, IT<:Integer}

Mesh of a quarter circle with a given number of elements per radius.

The parameter nperradius should be an even number; if that isn't so is adjusted to by adding one.

Q4elliphole(
    xradius::T,
    yradius::T,
    L::T,
    H::T,
    nL::IT,
    nH::IT,
    nW::IT,
) where {T<:Number, IT<:Integer}

Mesh of one quarter of a rectangular plate with an elliptical hole.

xradius,yradius = radius of the ellipse, L,H= and dimensions of the plate, nL,nH= numbers of edges along the side of the plate; this also happens to be the number of edges along the circumference of the elliptical hole nW= number of edges along the remaining straight edge (from the hole in the direction of the length),

Q4extrudeL2(
    fens::FENodeSet,
    fes::FESetL2,
    nLayers::IT,
    extrusionh::F,
) where {F<:Function, IT<:Integer}

Extrude a mesh of linear segments into a mesh of quadrilaterals (Q4).

Q4quadrilateral(xyz::Matrix{T}, nL::IT, nW::IT) where {T<:Number, IT<:Integer}

Mesh of a general quadrilateral given by the location of the vertices.

Q4refine(fens::FENodeSet, fes::FESetQ4)

Refine a mesh of quadrilaterals by bisection.

Q4spheren(radius::T, nperradius::IT) where {T<:Number, IT<:Integer}

Generate mesh of a spherical surface (1/8th of the sphere).

Q4toQ8(fens::FENodeSet, fes::FESetQ4)

Convert a mesh of quadrilateral Q4 to quadrilateral Q8.

Q8annulus(rin::T, rex::T, nr::IT, nc::IT, Angl::T) where {T<:Number, IT<:Integer}

Mesh of an annulus segment.

Mesh of an annulus segment, centered at the origin, with internal radius rin, and external radiusrex, and development angle Angl. Divided into elements:nr,nc` in the radial and circumferential direction respectively.

Q8block(Length::T, Width::T, nL::IT, nW::IT) where {T<:Number, IT<:Integer}

Mesh of a rectangle of Q8 elements.

Q8blockx(xs::Vector{T}, ys::Vector{T}) where {T<:Number, IT<:Integer}

Graded mesh of a 2-D block of Q8 finite elements.

Q9blockx(xs::Vector{T}, ys::Vector{T}) where {T<:Number}

Create a block of the quadratic Lagrangean Q9 nine-node quadrilaterals.

T10block(
    Length::T,
    Width::T,
    Height::T,
    nL::IT,
    nW::IT,
    nH::IT;
    orientation::Symbol = :a,
) where {T<:Number, IT<:Integer}

Generate a tetrahedral mesh of T10 elements of a rectangular block.

T10blockx(xs::VecOrMat{T}, ys::VecOrMat{T}, zs::VecOrMat{T}, orientation::Symbol = :a) where {T<:Number}

Generate a graded 10-node tetrahedral mesh of a 3D block.

10-node tetrahedra in a regular arrangement, with non-uniform given spacing between the nodes, with a given orientation of the diagonals.

The mesh is produced by splitting each logical rectangular cell into six tetrahedra.

T10cylindern(R::T, L::T, nR::IT, nL::IT; orientation = :b) where {T<:Number, IT<:Integer}

Ten-node tetrahedron mesh of solid cylinder with given number of edges per radius.

The axis of the cylinder is along the Z axis.

Even though the orientation is controllable, for some orientations the mesh is highly distorted (:a, :ca, :cb). So a decent mesh can only be expected for the orientation :b (default).

T10layeredplatex(
    xs::VecOrMat{T},
    ys::VecOrMat{T},
    ts::VecOrMat{T},
    nts::VecOrMat{IT},
    orientation::Symbol = :a,
) where {T<:Number, IT<:Integer}

T10 mesh for a layered block (composite plate) with specified in plane coordinates.

xs,ys =Locations of the individual planes of nodes. ts= Array of layer thicknesses, nts= array of numbers of elements per layer

The finite elements of each layer are labeled with the layer number, starting from 1 at the bottom.

T10quartercyln(R::T, L::T, nR::IT, nL::IT; orientation = :b) where {T<:Number, IT<:Integer}

Ten-node tetrahedron mesh of one quarter of solid cylinder with given number of edges per radius.

See: T4quartercyln

T10refine(fens::FENodeSet, fes::FESetT10)

Refine the mesh of quadratic tetrahedra.

Each tetrahedron is converted to eight tetrahedra (each face is quadri-sected).

T10toT4(fens::FENodeSet,  fes::FESetT4)

Convert a mesh of tetrahedra of type T10 (quadratic 10-node) to tetrahedra T4.

If desired, the array of nodes may be compacted so that unconnected nodes are deleted.

T4block(
    Length::T,
    Width::T,
    Height::T,
    nL::IT,
    nW::IT,
    nH::IT,
    orientation::Symbol = :a,
) where {T<:Number, IT<:Integer}

Generate a tetrahedral mesh of the 3D block.

Four-node tetrahedra in a regular arrangement, with uniform spacing between the nodes, with a given orientation of the diagonals.

The mesh is produced by splitting each logical rectangular cell into five or six tetrahedra, depending on the orientation: orientation = :a, :b generates 6 tetrahedra per cell. :ca, :cb generates 5 tetrahedra per cell.

Range =<0, Length> x <0, Width> x <0, Height>. Divided into elements: nL x nW x nH.

T4blockx(xs::VecOrMat{T}, ys::VecOrMat{T}, zs::VecOrMat{T}, orientation::Symbol = :a) where {T<:Number}

Generate a graded tetrahedral mesh of a 3D block.

Four-node tetrahedra in a regular arrangement, with non-uniform given spacing between the nodes, with a given orientation of the diagonals.

The mesh is produced by splitting each logical rectangular cell into five or six tetrahedra: refer to T4block.

T4cylindern(R::T, L::T, nR::IT, nL::IT; orientation = :b) where {T<:Number, IT<:Integer}

Four-node tetrahedron mesh of solid cylinder with given number of edges per radius.

The axis of the cylinder is along the Z axis.

Even though the orientation is controllable, for some orientations the mesh is highly distorted (:a, :ca, :cb). So a decent mesh can only be expected for the orientation :b (default).

T4extrudeT3(
    fens::FENodeSet,
    fes::FESetT3,
    nLayers::IT,
    extrusionh::F,
) where {F<:Function, IT<:Integer}

Extrude a mesh of triangles into a mesh of tetrahedra (T4).

T4meshedges(t::Array{IT,2}) where {IT<:Integer}

Compute all the edges of the 4-node triangulation.

T4quartercyln(R::T, L::T, nR::IT, nL::IT; orientation = :b) where {T<:Number, IT<:Integer}

Four-node tetrahedron mesh of one quarter of solid cylinder with given number of edges per radius.

The axis of the cylinder is along the Z axis. The mesh may be mirrored to create half a cylinder or a full cylinder.

Even though the orientation is controllable, for some orientations the mesh is highly distorted (:a, :ca, :cb). So a decent mesh can only be expected for the orientation :b (default).

T4refine(fens::FENodeSet, fes::FESetT4)

Refine a mesh of 4-node tetrahedra by octasection.

T4refine20(fens::FENodeSet, fes::FESetT4)

Refine a tetrahedral four-node mesh into another four-node tetrahedral mesh, with each original tetrahedron being subdivided into 20 new tetrahedra.

Each vertex is given one hexahedron. The scheme generates 15 nodes per tetrahedron when creating the hexahedra, one for each edge, one for each face, and one for the interior.

T4toT10(fens::FENodeSet,  fes::FESetT4)

Convert a mesh of tetrahedra of type T4 (four-node) to tetrahedra T10.

T4voximg(
    img::Array{DataT,3},
    voxdims::T,
    voxval::Array{DataT,1},
) where {DataT<:Number, T}

Generate a tetrahedral mesh from three-dimensional image.

tetv(X)

Compute the volume of a tetrahedron.

X = [0  4  3
9  2  4
6  1  7
0  1  5] # for these points the volume is 10.0
tetv(X)

tetv(X)

Compute the volume of a tetrahedron.

X = [0  4  3
9  2  4
6  1  7
0  1  5] # for these points the volume is 10.0
tetv(X)

tetv(X)

Compute the volume of a tetrahedron.

X = Float64[0  4  3
9  2  4
6  1  7
0  1  5] # for these points the volume is 10.0
tetv(X, 1, 2, 3, 4)

tetv1times6(v, i1, i2, i3, i4)

Compute 6 times the volume of the tetrahedron.

H20block(Length::T, Width::T, Height::T, nL::IT, nW::IT, nH::IT) where {T<:Number, IT<:Integer}

Create mesh of a 3-D block of H20 finite elements.

H20blockx(xs::Vector{T}, ys::Vector{T}, zs::Vector{T}) where {T<:Number}

Graded mesh of a 3-D block of H20 finite elements.

H27block(Length::T, Width::T, Height::T, nL::IT, nW::IT, nH::IT) where {T<:Number, IT<:Integer}

Create mesh of a 3-D block of H27 finite elements.

H27blockx(xs::Vector{T}, ys::Vector{T}, zs::Vector{T}) where {T<:Number}

Graded mesh of a 3-D block of H27 finite elements.

H8block(Length::T, Width::T, Height::T, nL::IT, nW::IT, nH::IT) where {T<:Number, IT<:Integer}

Make a mesh of a 3D block consisting of eight node hexahedra.

Length, Width, Height= dimensions of the mesh in Cartesian coordinate axes, smallest coordinate in all three directions is 0 (origin) nL, nW, nH=number of elements in the three directions

H8blockx(xs::Vector{T}, ys::Vector{T}, zs::Vector{T}) where {T<:Number}

Graded mesh of a 3-D block of H8 finite elements.

H8cylindern(Radius::T, Length::T, nperradius::IT, nL::IT) where {T<:Number, IT<:Integer}

H8 mesh of a solid cylinder with given number of edges per radius (nperradius) and per length (nL).

H8elliphole(
    xradius::T,
    yradius::T,
    L::T,
    H::T,
    T::T,
    nL::IT,
    nH::IT,
    nW::IT,
    nT::IT,
) where {T<:Number, IT<:Integer}

Mesh of one quarter of a rectangular plate with an elliptical hole.

xradius,yradius = radii of the ellipse, L,H = dimensions of the plate, Th = thickness of the plate nL,nH= numbers of edges along the side of the plate; this is also the number of edges along the circumference of the elliptical hole nW = number of edges along the remaining straight edge (from the hole in the radial direction),

H8extrudeQ4(
    fens::FENodeSet,
    fes::FESetQ4,
    nLayers::IT,
    extrusionh::F,
) where {F<:Function, IT<:Integer}

Extrude a mesh of quadrilaterals into a mesh of hexahedra (H8).

H8hexahedron(xyz::Matrix{T}, nL::IT, nW::IT, nH::IT; blockfun = nothing) where {T<:Number, IT<:Integer}

Mesh of a general hexahedron given by the location of the vertices.

xyz = One vertex location per row; Either two rows (for a rectangular block given by the its corners), or eight rows (general hexahedron). nL, nW, nH = number of elements in each direction blockfun = Optional argument: function of the block-generating mesh function (having the signature of the function H8block()).

H8layeredplatex(xs::Vector{T}, ys::Vector{T}, ts::Vector{T}, nts::Vector{IT}) where {T<:Number, IT<:Integer}

H8 mesh for a layered block (composite plate) with specified in plane coordinates.

xs,ys =Locations of the individual planes of nodes. ts= Array of layer thicknesses, nts= array of numbers of elements per layer

The finite elements of each layer are labeled with the layer number, starting from 1.

H8refine(fens::FENodeSet,  fes::FESetH8)

Refine a mesh of H8 hexahedrals by octasection.

H8sphere(radius::T, nrefine::IT) where {T<:Number, IT<:Integer}

Create a mesh of 1/8 of the sphere of "radius". The  mesh will consist of
four hexahedral elements if "nrefine==0",  or more if "nrefine>0".
"nrefine" is the number of bisections applied  to refine the mesh.

H8spheren(radius::T, nperradius::IT) where {T<:Number, IT<:Integer}

Create a solid mesh of 1/8 of sphere.

Create a solid mesh of 1/8 of the sphere of radius, with nperradius elements per radius.

H8toH20(fens::FENodeSet,  fes::FESetH8)

Convert a mesh of hexahedra H8 to hexahedra H20.

H8toH27(fens::FENodeSet,  fes::FESetH8)

Convert a mesh of hexahedra H8 to hexahedra H27.

H8voximg(
    img::Array{DataT,3},
    voxdims::Vector{T},
    voxval::Array{DataT,1},
) where {DataT<:Number}

Generate a hexahedral mesh from three-dimensional image.

T4toH8(fens::FENodeSet, fes::FESetT4)

Convert a tetrahedral four-node mesh into eight-node hexahedra.

Each vertex is given one hexahedron. The scheme generates 15 nodes per tetrahedron when creating the hexahedra, one for each edge, one for each face, and one for the interior.

compactnodes(fens::FENodeSet, connected::Vector{Bool})

Compact the finite element node set by deleting unconnected nodes.

fens = array of finite element nodes connected = The array element connected[j] is either 0 (when j is an unconnected node), or a positive number (when node j is connected to other nodes by at least one finite element)

Output

fens = new set of finite element nodes new_numbering= array which tells where in the new fens array the connected nodes are (or 0 when the node was unconnected). For instance, node 5 was connected, and in the new array it is the third node: then new_numbering[5] is 3.

Examples

Let us say there are nodes not connected to any finite element that you would like to remove from the mesh: here is how that would be accomplished.

connected = findunconnnodes(fens, fes);
fens, new_numbering = compactnodes(fens, connected);
fes = renumberconn!(fes, new_numbering);

Finally, check that the mesh is valid:

validate_mesh(fens, fes);

distortblock(
    B::F,
    Length::T,
    Width::T,
    nL::IT,
    nW::IT,
    xdispmul::T,
    ydispmul::T,
) where {F<:Function, T<:Number, IT<:Integer}

Distort a block mesh by shifting around the nodes. The goal is to distort the horizontal and vertical mesh lines into slanted lines. This is useful when testing finite elements where special directions must be avoided.

distortblock(ofens::FENodeSet{T}, xdispmul::T, ydispmul::T) where {T<:Number}

Distort a block mesh by shifting around the nodes. The goal is to distort the horizontal and vertical mesh lines into slanted lines.

element_coloring!(element_colors, unique_colors, color_counts, fes, n2e, ellist)

Find coloring of the elements such that no two elements of the same color share a node.

Compute element coloring, supplying the current element colors and the list of elements to be assigned colors.

element_coloring(fes, n2e, ellist::Vector{IT} = IT[]) where {IT<:Integer}

Find coloring of the elements such that no two elements of the same color share a node.

Returns the colors of the elements (color here means an integer), and a list of the unique colors (numbers).

ellist = list of elements to be assigned colors; other element colors may be looked at

fusenodes(fens1::FENodeSet{T}, fens2::FENodeSet{T}, tolerance::T) where {T<:Number}

Fuse together nodes from two node sets.

Fuse two node sets. If necessary, by gluing together nodes located within tolerance of each other. The two node sets, fens1 and fens2, are fused together by merging the nodes that fall within a box of size tolerance. The merged node set, fens, and the new indexes of the nodes in the set fens1 are returned.

The set fens2 will be included unchanged, in the same order, in the node set fens. The indexes of the node set fens1 will have changed.

Example

After the call to this function we have k=new_indexes_of_fens1_nodes[j] is the node in the node set fens which used to be node j in node set fens1. The finite element set connectivity that used to refer to fens1 needs to be updated to refer to the same nodes in the set fens as updateconn!(fes, new_indexes_of_fens1_nodes);

interior2boundary(interiorconn::Array{IT,2}, extractb::Array{IT,2}) where {IT<:Integer}

Extract the boundary connectivity from the connectivity of the interior.

extractb = array that defines in which order the bounding faces are traversed. For example, for tetrahedra this array is extractb = [1 3 2; 1 2 4; 2 3 4; 1 4 3]

mergemeshes(
    fens1::FENodeSet{T},
    fes1::T1,
    fens2::FENodeSet{T},
    fes2::T2,
    tolerance::T,
) where {T, T1<:AbstractFESet, T2<:AbstractFESet}

Merge together two meshes.

Merge two meshes together by gluing together nodes within tolerance. The two meshes, fens1, fes1, and fens2, fes2, are glued together by merging the nodes that fall within a box of size tolerance. If tolerance is set to zero, no merging of nodes is performed; the two meshes are simply concatenated together.

The merged node set, fens, and the two finite element sets with renumbered connectivities are returned.

Important notes: On entry into this function the connectivity of fes1 point into fens1 and the connectivity of fes2 point into fens2. After this function returns the connectivity of both fes1 and fes2 point into fens. The order of the nodes of the node set fens1 in the resulting set fens will have changed, whereas the order of the nodes of the node set fens2 is are guaranteed to be the same. Therefore, the connectivity of fes2 will in fact remain the same.

mergenmeshes(meshes::Array{Tuple{FENodeSet{T},AbstractFESet}}, tolerance::T) where {T<:Number}

Merge several meshes together.

The meshes are glued together by merging the nodes that fall within a box of size tolerance. If tolerance is set to zero, no merging of nodes is performed; the nodes from the meshes are simply concatenated together.

Output

The merged node set, fens, and an array of finite element sets with renumbered connectivities are returned.

mergenodes(
    fens::FENodeSet{T},
    fes::AbstractFESet,
    tolerance::T,
    candidates::AbstractVector{IT},
) where {T<:Number, IT<:Integer}

Merge together nodes of a single node set.

Similar to mergenodes(fens, fes, tolerance), but only the candidate nodes are considered for merging. This can potentially speed up the operation by orders of magnitude.

mergenodes(fens::FENodeSet{T}, fes::AbstractFESet, tolerance::T) where {T<:Number}

Merge together nodes of a single node set.

Merge by gluing together nodes from a single node set located within tolerance of each other. The nodes are glued together by merging the nodes that fall within a box of size tolerance. The merged node set, fens, and the finite element set, fes, with renumbered connectivities are returned.

Warning: This tends to be an expensive operation!

meshboundary(fes::T) where {T<:AbstractFESet}

Compute the boundary of a mesh defined by the given finite element set.

Arguments

• fes::T: The finite element set representing the mesh.

Returns

The boundary of the mesh.

Extract the finite elements of manifold dimension (n-1) from the supplied finite element set of manifold dimension (n).

meshsmoothing(fens::FENodeSet, fes::T; options...) where {T<:AbstractFESet}

General smoothing of meshes.

Keyword options

method = :taubin (default) or :laplace fixedv = Boolean array, one entry per vertex: is the vertex immovable (true) or movable (false) npass = number of passes (default 2)

Return

The modified node set.

mirrormesh(
    fens::FENodeSet,
    fes::ET,
    Normal::Vector{T},
    Point::Vector{T};
    kwargs...,
) where {ET<:AbstractFESet, T<:Number}

Mirror a mesh in a plane given by its normal and one point.

Keyword arguments

• renumb = node renumbering function for the mirrored element

Warning: The code to relies on the numbering of the finite elements: to reverse the orientation of the mirrored finite elements, the connectivity is listed in reverse order. If the mirrored finite elements do not follow this rule (for instance hexahedra or quadrilaterals), their areas/volumes will come out negative. In such a case the renumbering function of the connectivity needs to be supplied.

For instance: H8 elements require the renumbering function to be supplied as

renumb = (c) -> c[[1, 4, 3, 2, 5, 8, 7, 6]]

nodepartitioning(fens::FENodeSet, nincluded::Vector{Bool}, npartitions)

Compute the inertial-cut partitioning of the nodes.

nincluded = Boolean array: is the node to be included in the partitioning or not? npartitions = number of partitions, but note that the actual number of partitions is going to be a power of two.

The partitioning can be visualized for instance as:

partitioning = nodepartitioning(fens, npartitions)
partitionnumbers = unique(partitioning)
for gp = partitionnumbers
  groupnodes = findall(k -> k == gp, partitioning)
  File =  "partition-nodes-$(gp).vtk"
  vtkexportmesh(File, fens, FESetP1(reshape(groupnodes, length(groupnodes), 1)))
end
File =  "partition-mesh.vtk"
vtkexportmesh(File, fens, fes)
@async run(`"paraview.exe" $File`)

nodepartitioning(fens::FENodeSet, npartitions)

Compute the inertial-cut partitioning of the nodes.

npartitions = number of partitions, but note that the actual number of partitions will be a power of two.

In this variant all the nodes are to be included in the partitioning.

The partitioning can be visualized for instance as:

partitioning = nodepartitioning(fens, npartitions)
partitionnumbers = unique(partitioning)
for gp = partitionnumbers
  groupnodes = findall(k -> k == gp, partitioning)
  File =  "partition-nodes-Dollar(gp).vtk"
  vtkexportmesh(File, fens, FESetP1(reshape(groupnodes, length(groupnodes), 1)))
end
File =  "partition-mesh.vtk"
vtkexportmesh(File, fens, fes)
@async run(`"paraview.exe" DollarFile`)

nodepartitioning(fens::FENodeSet, fesarr, npartitions::Vector{Int})

Compute the inertial-cut partitioning of the nodes.

fesarr = array of finite element sets that represent regions npartitions = array of the number of partitions in each region. However note that the actual number of partitions will be a power of two.

The partitioning itself is carried out by nodepartitioning() with a list of nodes to be included in the partitioning. For each region I the nodes included in the partitioning are those connected to the elements of that region, but not to elements that belong to any of the previous regions, 1…I-1.

outer_surface_of_solid(fens, bdry_fes)

Find the finite elements that form the outer boundary surface.

!!! note:

The code will currently not work correctly for 2D axially symmetric geometries.

Return

Set of boundary finite elements that enclose the solid. Now cavities are included.

pointpartitioning(xyz, npartitions::Int = 2)

Compute the inertial-cut partitioning of a set of points.

renumberconn!(fes::AbstractFESet, new_numbering::AbstractVector{IT}) where {IT<:Integer}

Renumber the nodes in the connectivity of the finite elements based on a new numbering for the nodes.

fes =finite element set new_numbering = new serial numbers for the nodes. The connectivity should be changed as conn[j] –> new_numbering[conn[j]]

Let us say there are nodes not connected to any finite element that we would like to remove from the mesh: here is how that would be accomplished.

connected = findunconnnodes(fens, fes);
fens, new_numbering = compactnodes(fens, connected);
fes = renumberconn!(fes, new_numbering);

Finally, check that the mesh is valid:

validate_mesh(fens, fes);

reordermesh(fens, fes, ordering)

Reorder mesh (reshuffle nodes, renumber connectivities correspondingly).

Arguments

• fens: The set of mesh nodes.
• fes: The set of elements.
• ordering: The desired ordering of the nodes and elements.

Returns

The reordered mesh nodes and elements.

The ordering may come from Reverse Cuthill-McKey (package SymRCM).

vertexneighbors(conn::Matrix{IT}, nvertices::IT) where {IT<:Integer}

Find the node neighbors in the mesh.

Return an array of integer vectors, element I holds an array of numbers of nodes which are connected to node I (including node I).

vsmoothing(v::Matrix{T}, t::Matrix{IT}; kwargs...) where {T<:Number, IT<:Integer}

Internal routine for mesh smoothing.

Keyword options: method = :taubin (default) or :laplace fixedv = Boolean array, one entry per vertex: is the vertex immovable (true) or movable (false) npass = number of passes (default 2)

addhyperface!(container,hyperface,newn)

Add a hyper face to the container.

The new node is stored in hyper face data in the container and can be retrieved later.

findhyperface!(container,hyperface)

Find a hyper face in the container.

If the container holds the indicated hyper face, it is returned together with the stored new node number.

gradedspace(start::T, stop::T, N::Int)  where {T<:Number}

Generate quadratic space.

Generate a quadratic sequence of N numbers between start and stop. This sequence corresponds to separation of adjacent numbers that increases linearly from start to finish.

Example

julia> gradedspace(2.0, 3.0, 5)
5-element Array{Float64,1}:
 2.0
 2.0625
 2.25
 2.5625
 3.0

linearspace(start::T, stop::T, N::Int)  where {T<:Number}

Generate linear space.

Generate a linear sequence of N numbers between start and stop (i. e. sequence of number with uniform intervals inbetween).

Example

julia> linearspace(2.0, 3.0, 5)
2.0:0.25:3.0

logspace(start::T, stop::T, N::Int)  where {T<:Number}

Generate logarithmic space.

Generate a logarithmic sequence of N numbers between start and stop (i. e. sequence of number with uniform intervals inbetween in the log space).

Example

julia> logspace(2.0, 3.0, 5)                                                             
5-element Array{Float64,1}:                                                              
  100.0
  177.82794100389228   
  316.2277660168379      
  562.341325190349  
 1000.0    

makecontainer()

Make hyper face container.

This is a dictionary of hyper faces, indexed with an anchor node. The anchor node is the smallest node number within the connectivity of the hyper face.

ontosphere(xyz::Matrix{T}, radius::T) where {T}

Project nodes onto a sphere of given radius.

Arguments

• xyz::Matrix{T}: A matrix of shape (3, N) representing the coordinates of N points.
• radius::T: The radius of the sphere.

Returns

A matrix of shape (3, N) representing the coordinates of points on the sphere.

import_ABAQUS(filename)

Import Abaqus mesh (.inp file).

Limitations:

1. Only the *NODE and *ELEMENT sections are read
2. Only 4-node and 10-node tetrahedra, 8-node or 20-node hexahedra, 4-node quadrilaterals, 3-node triangles are handled.

Output

Data dictionary, with keys

• "fens" = finite element nodes.
• "fesets" = array of finite element sets.
• "nsets" = dictionary of "node sets", the keys are the names of the sets.

import_H5MESH(meshfile)

Import mesh in the H5MESH format (.h5mesh file).

Output

Data dictionary, with keys

• "fens" = finite element nodes.
• "fesets" = array of finite element sets.

import_MESH(filename)

Import mesh in the MESH format (.mesh, .xyz, .conn triplet of files).

Output

Data dictionary, with keys

• "fens" = finite element nodes.
• "fesets" = array of finite element sets.

import_NASTRAN(filename)

Import tetrahedral (4- and 10-node) NASTRAN mesh (.nas file).

Limitations:

1. only the GRID, CTETRA, and PSOLID sections are read.
2. Only 4-node and 10-node tetrahedra are handled.
3. The file should be free-form (data separated by commas).

Some fixed-format files can also be processed (large-field, but not small-field).

Output

Data dictionary, with keys

• "fens": set of finite element nodes,
• "fesets": array of finite element sets,
• "property_ids": array of property ids (integers) associated with the finite element sets.

vtkexportmesh(theFile::String, Connectivity, Points, Cell_type;
    vectors=nothing, scalars=nothing)

Export mesh to a VTK 1.0 file as an unstructured grid.

Arguments:

• theFile = file name,
• Connectivity = array of connectivities, one row per element,
• Points = array of node coordinates, one row per node,
• Cell_type = type of the cell, refer to the predefined constants P1, L2, ..., H20, ...
• scalars = array of tuples, (name, data)
• vectors = array of tuples, (name, data)

For the scalars: If data is a vector, that data is exported as a single field. On the other hand, if it is an 2d array, each column is exported as a separate field.

vtkexportmesh(theFile::String, fens::FENodeSet, fes::T; opts...) where {T<:AbstractFESet}

Export mesh to a VTK 1.0 file as an unstructured grid.

Arguments:

• theFile = file name,
• fens = finite element node set,
• fes = finite element set,
• opts = keyword argument list, where
  • scalars = array of tuples, (name, data)
  • vectors = array of tuples, (name, data)

For the scalars: If data is a vector, that data is exported as a single field. On the other hand, if it is an 2d array, each column is exported as a separate field.

vtkexportvectors(theFile::String, Points, vectors)

Export vector data to a VTK 1.0 file.

Arguments:

• theFile = file name,
• Points = array of collection of coordinates (tuples or vectors),
• vectors = array of tuples, (name, data), where name is a string, and data is array of collection of coordinates (tuples or vectors).

Example

Points = [(1.0, 3.0), (0.0, -1.0)]
vectors = [("v", [(-1.0, -2.0), (1.0, 1.0)])]
vtkexportvectors("theFile.VTK", Points, vectors)

will produce file with

# vtk DataFile Version 1.0
FinEtools Export
ASCII

DATASET UNSTRUCTURED_GRID
POINTS 2 double
1.0 3.0 0.0
0.0 -1.0 0.0


POINT_DATA 2
VECTORS v double
-1.0 -2.0 0.0
1.0 1.0 0.0

close(self::AbaqusExporter)

Close the stream opened by the exporter.

ASSEMBLY(self::AbaqusExporter, NAME::AbstractString)

Write out the *ASSEMBLY option.

BOUNDARY(self::AbaqusExporter, NSET::AbstractString, dof::Integer,  fixed_value)

Write out the *BOUNDARY option.

• NSET = name of a node set,
• is_fixed= array of Boolean flags (true means fixed, or prescribed), one row per node,
• fixed_value=array of displacements to which the corresponding displacement components is fixed

BOUNDARY(self::AbaqusExporter, NSET::AbstractString, dof::Integer)

Write out the *BOUNDARY option to fix displacements at zero for a node set.

Invoke at Level: Model, Step

• NSET= node set,
• dof=Degree of freedom, 1, 2, 3

BOUNDARY(self::AbaqusExporter, NSET::AbstractString, dof::Integer,
  value::F) where {F}

Write out the *BOUNDARY option to fix displacements at nonzero value for a node set.

• NSET= node set,
• dof=Degree of freedom, 1, 2, 3
• typ = DISPLACEMENT

BOUNDARY(self::AbaqusExporter, meshornset, is_fixed::AbstractArray{B,2}, fixed_value::AbstractArray{F,2}) where {B, F}

Write out the *BOUNDARY option.

• meshornset = name of a mesh or a node set,
• is_fixed= array of Boolean flags (true means fixed, or prescribed), one row per node, as many columns as there are degrees of freedom per node,
• fixed_value=array of displacements to which the corresponding displacement components is fixed, as many columns as there are degrees of freedom per node

BOUNDARY(self::AbaqusExporter, mesh, nodes, is_fixed::AbstractArray{B,2}, fixed_value::AbstractArray{F,2}) where {B, F}

Write out the *BOUNDARY option.

The boundary condition is applied to the nodes specified in the array nodes, in the mesh (or node set) meshornset.

meshornset = mesh or node set (default is empty) nodes = array of node numbers, the node numbers are attached to the mesh label, is_fixed= array of Boolean flags (true means fixed, or prescribed), one row per node, fixed_value=array of displacements to which the corresponding displacement components is fixed

Example

BOUNDARY(AE, "ASSEM1.INSTNC1", 1:4, fill(true, 4, 1), reshape([uy(fens.xyz[i, :]...) for i in 1:4], 4, 1))

CLOAD(self::AbaqusExporter, NSET::AbstractString, dof::Integer,
  magnitude::F) where {F}

Write out the *CLOAD option.

NSET=Node set dof= 1, 2, 3, magnitude= signed multiplier

CLOAD(self::AbaqusExporter, nodenumber::Integer, dof::Integer,
  magnitude::F) where {F}

Write out the *CLOAD option.

nodenumber=Number of node dof= 1, 2, 3, magnitude= signed multiplier

COMMENT(self::AbaqusExporter, Text::AbstractString)

Write out the ** comment option.

DENSITY(self::AbaqusExporter, rho)

Write out the *DENSITY option.

DLOAD(self::AbaqusExporter, ELSET::AbstractString,
  traction::AbstractVector{F}) where {F}

Write out the *DLOAD option.

ELASTIC(self::AbaqusExporter, E::F, nu::F) where {F}

Write out the *ELASTIC,TYPE=ISOTROPIC option.

ELASTIC(self::AbaqusExporter, E1::F, E2::F, E3::F, nu12::F, nu13::F, nu23::F,
  G12::F, G13::F, G23::F) where {F}

Write out the *ELASTIC,TYPE=ENGINEERING CONSTANTS option.

ELEMENT(self::AbaqusExporter, TYPE::AbstractString, ELSET::AbstractString,
  start::Integer, conn::AbstractArray{T, 2}) where {T<:Integer}

Write out the *ELEMENT option.

TYPE= element type code, ELSET= element set to which the elements belong, start= start the element numbering at this integer, conn= connectivity array for the elements, one row per element

ELSET_ELSET(self::AbaqusExporter, ELSET::AbstractString, n::AbstractArray{T, 1}) where {T<:Integer}

Write out the *ELSET option.

ELSET = name of the set, n = array of the node numbers

EL_PRINT(self::AbaqusExporter, ELSET::AbstractString, KEYS::AbstractString)

Write out the *EL PRINT option.

END_ASSEMBLY(self::AbaqusExporter)

Write out the *END ASSEMBLY option.

END_INSTANCE(self::AbaqusExporter)

Write out the *END INSTANCE option.

END_PART(self::AbaqusExporter)

Write out the *END PART option.

END_STEP(self::AbaqusExporter)

Write out the *END STEP option.

ENERGY_PRINT(self::AbaqusExporter)

Write out the *ENERGY PRINT option.

EXPANSION(self::AbaqusExporter, CTE::F) where {F}

Write out the *EXPANSION option.

HEADING(self::AbaqusExporter, Text::AbstractString)

Write out the *HEADING option.

INSTANCE(self::AbaqusExporter, NAME::AbstractString, PART::AbstractString)

Write out the *INSTANCE option.

MATERIAL(self::AbaqusExporter, MATERIAL::AbstractString)

Write out the *MATERIAL option.

NODE(self::AbaqusExporter, xyz::AbstractArray{T, 2}) where {T}

Write out the *NODE option.

xyz=array of node coordinates

NODE_PRINT(self::AbaqusExporter, NSET::AbstractString)

Write out the *NODE PRINT option.

NSET_NSET(self::AbaqusExporter, NSET::AbstractString,
  n::AbstractVector{T}) where {T<:Integer}

Write out the *NSET option.

NSET = name of the set, n = array of the node numbers

ORIENTATION(self::AbaqusExporter, ORIENTATION::AbstractString,
  a::AbstractArray{T,1}, b::AbstractArray{T,1})

Write out the *ORIENTATION option.

Invoke at level: Part, Part instance, Assembly

PART(self::AbaqusExporter, NAME::AbstractString)

Write out the *PART option.

SECTION_CONTROLS(self::AbaqusExporter, NAME::AbstractString,
  OPTIONAL::AbstractString)

Write out the *SECTION CONTROLS option.

OPTIONAL = string, for instance HOURGLASS=ENHANCED

SOLID_SECTION(self::AbaqusExporter, MATERIAL::AbstractString,
  ORIENTATION::AbstractString, ELSET::AbstractString)

Write out the *SOLID SECTION option.

Level: Part, Part instance

SOLID_SECTION(self::AbaqusExporter, MATERIAL::AbstractString,
  ORIENTATION::AbstractString, ELSET::AbstractString,
  CONTROLS::AbstractString)

Write out the *SOLID SECTION option.

Level: Part, Part instance

SOLID_SECTION(self::AbaqusExporter, MATERIAL::AbstractString,
  ORIENTATION::AbstractString, ELSET::AbstractString)

Write out the *SOLID SECTION option.

Level: Part, Part instance

STEP_FREQUENCY(self::AbaqusExporter, Nmodes::Integer)

Write out the *STEP,FREQUENCY option.

STEP_PERTURBATION_BUCKLE(self::AbaqusExporter, neigv::Integer)

Write out the *STEP,PERTURBATION option for linear buckling analysis.

STEP_PERTURBATION_STATIC(self::AbaqusExporter)

Write out the *STEP,PERTURBATION option for linear static analysis.

SURFACE_SECTION(self::AbaqusExporter, ELSET::AbstractString)

Write out the *SURFACE SECTION option.

TEMPERATURE(self::AbaqusExporter, nlist::AbstractArray{I, 1},
  tlist::AbstractArray{F, 1}) where {I, F}

Write out the *TEMPERATURE option.

close(self::NASTRANExporter)

Close the stream opened by the exporter.

BEGIN_BULK(self::NASTRANExporter)

Terminate the Case Control section by starting the bulk section.

CEND(self::NASTRANExporter)

Terminate the Executive Control section.

CTETRA(self::NASTRANExporter, eid::Int, pid::Int, conn::Vector{Int})

Write a statement for a single tetrahedron element.

ENDDATA(self::NASTRANExporter)

Terminate the bulk section.

GRID(self::NASTRANExporter, n::Int, xyz)

Write a grid-point statement.

MAT1(
    self::NASTRANExporter,
    mid::Int,
    E::T,
    nu::T,
    rho::T = 0.0,
    A::T = 0.0,
    TREF::T = 0.0,
    GE::T = 0.0,
) where {T}

Write a statement for an isotropic elastic material.

MAT1(
    self::NASTRANExporter,
    mid::Int,
    E::T,
    G::T,
    nu::T,
    rho::T,
    A::T,
    TREF::T,
    GE::T,
) where {T}

Write a statement for an isotropic elastic material.

PSOLID(self::NASTRANExporter, pid::Int, mid::Int)

Write solid-property statement.

close(self::STLExporter)

Close the stream opened by the exporter.

endsolid(self::STLExporter, name::AbstractString = "thesolid")

Write statement to end the solid.

facet(self::STLExporter, v1::Vector{T}, v2::Vector{T}, v3::Vector{T}) where {T}

Write a single facet.

solid(self::STLExporter, name::AbstractString = "thesolid")

Write a statement to begin the solid.

savecsv(name::String; kwargs...)

Save arrays as a CSV file.

Example:

savecsv("ab", a = rand(3), b = rand(3))

h2libexporttri(theFile::String, Connectivity, Points)

Write a file in the H2Lib format.

vtkwrite(theFile::String, Connectivity, Points, celltype; vectors=nothing, scalars=nothing)

Export mesh to VTK as an unstructured grid (binary format).

Arguments:

• theFile = file name,
• Connectivity = array of connectivities, one row per element,
• Points = array of node coordinates, one row per node,
• Cell_type = type of the cell, refer to the predefined constants WriteVTK.P1, WriteVTK.L2, ..., WriteVTK.H20`, ...
• scalars = array of tuples, (name, data)
• vectors = array of tuples, (name, data)

For the scalars: If data is a vector, that data is exported as a single field. On the other hand, if it is an 2d array, each column is exported as a separate field.

Return the vtk file.

vtkwrite(theFile::String, fens::FENodeSet, fes::T; opts...) where {T<:AbstractFESet}

Export mesh to VTK as an unstructured grid (binary file).

Arguments:

• theFile = file name,
• fens = finite element node set,
• fes = finite element set,
• opts = keyword argument list, where scalars = array of tuples, (name, data), vectors = array of tuples, (name, data)

For the scalars: If data is a vector, that data is exported as a single field. On the other hand, if it is an 2d array, each column is exported as a separate field.

vtkwritecollection(theFile::String, Connectivity, Points, celltype, times; vectors=nothing, scalars=nothing)

Write a collection of VTK files (.pvd file).

times: array of times

All the other arguments are the same as for vtkwrite. If scalars or vectors are supplied, they correspond to the times in the times array.

See the vtkwritecollection methods.

vtkwritecollection(theFile::String, fens::FENodeSet, fes::T, times; opts...) where {T<:AbstractFESet}

Write a collection of VTK files (.pvd file).

times: array of times

All the other arguments are the same as for vtkwrite. If scalars or vectors are supplied, they correspond to the times in the times array.

See the vtkwritecollection methods.

write_H5MESH(meshfile::String, fens::FENodeSet, fes::T) where {T<:AbstractFESet}

Write the mesh in the H5MESH format.

The mesh is stored in a HDF5 file.

associategeometry!(self::AbstractFEMM, geom::NodalField{GFT}) where {GFT}

Associate geometry field with the FEMM.

There may be operations that could benefit from pre-computations that involve a geometry field. If so, associating the geometry field gives the FEMM a chance to save on repeated computations.

Geometry field is normally passed into any routine that evaluates some forms (integrals) over the mesh. Whenever the geometry passed into a routine is not consistent with the one for which associategeometry!() was called before, associategeometry!() needs to be called with the new geometry field.

bilform_convection(
    self::FEMM,
    assembler::A,
    geom::NodalField{FT},
    u::NodalField{T},
    Q::NodalField{QT},
    rhof::DC
) where {FEMM<:AbstractFEMM, A<:AbstractSysmatAssembler, FT, T, QT, DC<:DataCache}

Compute the sparse matrix implied by the bilinear form of the "convection" type.

\[\int_{V} {w} \rho \mathbf{u} \cdot \nabla q \; \mathrm{d} V\]

Here $w$ is the scalar test function, $\mathbf{u}$ is the convective velocity, $q$ is the scalar trial function, $\rho$ is the mass density; $\rho$ is computed by rhof, which is a given function(data). Both test and trial functions are assumed to be from the same approximation space. rhof is represented with DataCache, and needs to return a scalar mass density.

The integral is with respect to the volume of the domain $V$ (i.e. a three dimensional integral).

Arguments

• self = finite element machine;
• assembler = assembler of the global matrix;
• geom = geometry field;
• u = convective velocity field;
• Q = nodal field to define the degree of freedom numbers;
• rhof= data cache, which is called to evaluate the coefficient $\rho$, given the location of the integration point, the Jacobian matrix, and the finite element label.

bilform_diffusion(
    self::FEMM,
    assembler::A,
    geom::NodalField{FT},
    u::NodalField{T},
    cf::DC
) where {FEMM<:AbstractFEMM, A<:AbstractSysmatAssembler, FT, T, DC<:DataCache}

Compute the sparse matrix implied by the bilinear form of the "diffusion" type.

\[\int_{V} \nabla w \cdot c \cdot \nabla u \; \mathrm{d} V\]

Here $\nabla w$ is the gradient of the scalar test function, $\nabla u$ is the gradient of the scalar trial function, $c$ is a square symmetric matrix of coefficients or a scalar; $c$ is computed by cf, which is a given function (data). Both test and trial functions are assumed to be from the same approximation space. cf is represented with DataCache, and needs to return a symmetric square matrix (to represent general anisotropic diffusion) or a scalar (to represent isotropic diffusion).

The coefficient matrix $c$ can be given in the so-called local material coordinates: coordinates that are attached to a material point and are determined by a local cartesian coordinates system (mcsys).

The integral is with respect to the volume of the domain $V$ (i.e. a three dimensional integral).

Arguments

• self = finite element machine;
• assembler = assembler of the global matrix;
• geom = geometry field;
• u = nodal field to define the degree of freedom numbers;
• cf= data cache, which is called to evaluate the coefficient $c$, given the location of the integration point, the Jacobian matrix, and the finite element label.

bilform_div_grad(
    self::FEMM,
    assembler::A,
    geom::NodalField{FT},
    u::NodalField{T},
    viscf::DC
) where {FEMM<:AbstractFEMM, A<:AbstractSysmatAssembler, FT, T, DC<:DataCache}

Compute the sparse matrix implied by the bilinear form of the "div grad" type.

\[\int_{V} \mu \nabla \mathbf{w}: \nabla\mathbf{u} \; \mathrm{d} V\]

Here $\mathbf{w}$ is the vector test function, $\mathbf{u}$ is the velocity, $\mu$ is the dynamic viscosity (or kinematic viscosity, depending on the formulation); $\mu$ is computed by viscf, which is a given function (data). Both test and trial functions are assumed to be from the same approximation space. viscf is represented with DataCache, and needs to return a scalar viscosity.

The integral is with respect to the volume of the domain $V$ (i.e. a three dimensional integral).

Arguments

• self = finite element machine;
• assembler = assembler of the global matrix;
• geom = geometry field;
• u = velocity field;
• viscf= data cache, which is called to evaluate the coefficient $\mu$, given the location of the integration point, the Jacobian matrix, and the finite element label.

bilform_dot(
    self::FEMM,
    assembler::A,
    geom::NodalField{FT},
    u::NodalField{T},
    cf::DC
) where {FEMM<:AbstractFEMM, A<:AbstractSysmatAssembler, FT, T, DC<:DataCache}

Compute the sparse matrix implied by the bilinear form of the "dot" type.

\[\int_{\Omega} \mathbf{w} \cdot \mathbf{c} \cdot \mathbf{u} \; \mathrm{d} \Omega\]

Here $\mathbf{w}$ is the test function, $\mathbf{u}$ is the trial function, $\mathbf{c}$ is a square matrix of coefficients; $\mathbf {c}$ is computed by cf, which is a given function (data). Both trial and test functions are assumed to be vectors(even if of length 1). cf is represented with DataCache, and needs to return a square matrix, with dimension equal to the number of degrees of freedom per node in the u field.

The integral domain $\Omega$ can be the volume of the domain $V$ (i.e. a three dimensional integral), or a surface $S$ (i.e. a two dimensional integral), or a line domain $L$ (i.e. a one dimensional integral).

Arguments

• self = finite element machine;
• assembler = assembler of the global object;
• geom = geometry field;
• u = nodal field to define the degree of freedom numbers;
• cf= data cache, which is called to evaluate the coefficient $c$, given the location of the integration point, the Jacobian matrix, and the finite element label.
• m = manifold dimension (default is 3).

bilform_lin_elastic(
    self::FEMM,
    assembler::A,
    geom::NodalField{FT},
    u::NodalField{T},
    cf::DC
) where {FEMM<:AbstractFEMM, A<:AbstractSysmatAssembler, FT, T, DC<:DataCache}

Compute the sparse matrix implied by the bilinear form of the "linearized elasticity" type.

\[\int_{V} (B \mathbf{w})^T C B \mathbf{u} \; \mathrm{d} V\]

Here $\mathbf{w}$ is the vector test function, $\mathbf{u}$ is the displacement (velocity), $C$ is the elasticity (viscosity) matrix; $C$ is computed by cf, which is a given function(data). Both test and trial functions are assumed to be from the same approximation space. cf is represented with DataCache, and needs to return a matrix of the appropriate size.

The integral is with respect to the volume of the domain $V$ (i.e. a three dimensional integral).

Arguments

• self = finite element machine;
• assembler = assembler of the global matrix;
• geom = geometry field;
• u = velocity field;
• viscf= data cache, which is called to evaluate the coefficient $\mu$, given the location of the integration point, the Jacobian matrix, and the finite element label.

connectionmatrix(self::FEMM, nnodes) where {FEMM<:AbstractFEMM}

Compute the connection matrix.

The matrix has a nonzero in all the rows and columns which correspond to nodes connected by some finite element.

distribloads(
    self::FEMM,
    assembler::A,
    geom::NodalField{FT},
    P::NodalField{T},
    fi::ForceIntensity,
    m,
) where {FEMM<:AbstractFEMM, A<:AbstractSysvecAssembler, FT<:Number, T}

Compute distributed loads vector.

Arguments

• fi=force intensity object
• m= manifold dimension, 0= vertex (point), 1= curve, 2= surface, 3= volume. For body loads m is set to 3, for tractions on the surface it is set to 2, and so on.

The actual work is done by linform_dot().

dualconnectionmatrix(
    self::FEMM,
    fens::FENodeSet,
    minnodes = 1,
) where {FEMM<:AbstractFEMM}

Compute the dual connection matrix.

The matrix has a nonzero in all the rows and columns which correspond to elements connected by some finite element nodes.

• minnodes: minimum number of nodes that the elements needs to share in order to be neighbors (default 1)

elemfieldfromintegpoints(
    self::FEMM,
    geom::NodalField{GFT},
    u::NodalField{UFT},
    dT::NodalField{TFT},
    quantity::Symbol,
    component::AbstractVector{IT};
    context...,
) where {FEMM<:AbstractFEMM, GFT<:Number, UFT<:Number, TFT<:Number, IT<:Integer}

Construct elemental field from integration points.

Arguments

geom - reference geometry field u - displacement field dT - temperature difference field quantity - this is what you would assign to the 'quantity' argument of the material update!() method. component- component of the 'quantity' array: see the material update() method.

Output

• the new field that can be used to map values to colors and so on

ev_integrate(
        self::FEMM,
        geom::NodalField{FT},
        f::DC,
        initial::R,
        m,
) where {FEMM<:AbstractFEMM, FT<:Number, DC<:DataCache, R}

Compute the integral of a given function over a mesh domain.

\[\int_{\Omega} {f} \; \mathrm{d} \Omega\]

Here ${f}$ is a given function (data). The data ${f}$ is represented with DataCache.

Arguments

• self = finite element machine;
• geom = geometry field;
• f= data cache, which is called to evaluate the integrand based on the location, the Jacobian matrix, the finite element identifier, and the quadrature point;
• initial = initial value of the integral,
• m= manifold dimension, 0= vertex (point), 1= curve, 2= surface, 3= volume. For body loads m is set to 3, for tractions on the surface it is set to 2, and so on.

field_elem_to_nodal!(
    self::AbstractFEMM,
    geom::NodalField{FT},
    ef::EFL,
    nf::NFL;
    kind = :weighted_average,
) where {FT, T<:Number, EFL<:ElementalField{T}, NFL<:NodalField{T}}

Make a nodal field from an elemental field over the discrete manifold.

ef = ELEMENTAL field to supply the values nf = NODAL field to receive the values kind = default is :weighted_average; other options: :max

Returns nf.

field_nodal_to_elem!(
    self::AbstractFEMM,
    geom::NodalField{FT},
    nf::NFL,
    ef::EFL;
    kind = :weighted_average,
) where {FT<:Number, T, EFL<:ElementalField{T}, NFL<:NodalField{T}}

Make an elemental field from a nodal field over the discrete manifold.

nf = NODAL field to supply the values ef = ELEMENTAL field to receive the values kind = default is :weighted_average; other options: :max

Returns ef.

fieldfromintegpoints(
    self::FEMM,
    geom::NodalField{GFT},
    u::NodalField{UFT},
    dT::NodalField{TFT},
    quantity::Symbol,
    component::AbstractVector{IT};
    context...,
) where {FEMM<:AbstractFEMM, GFT<:Number, UFT<:Number, TFT<:Number, IT<:Integer}

Construct nodal field from integration points.

Arguments

• geom - reference geometry field
• u - displacement field
• dT - temperature difference field
• quantity - this is what you would assign to the 'quantity' argument of the material update!() method.
• component- component of the 'quantity' array: see the material update() method.

Keyword arguments

• nodevalmethod = :invdistance (the default) or :averaging;
• reportat = at which point should the element quantities be reported? This argument is interpreted inside the inspectintegpoints() method.

Output

• the new field that can be used to map values to colors and so on

finite_elements(self::FEMM) where {FEMM <: AbstractFEMM}

Retrieve the finite element set for this FEMM to work on.

innerproduct(
    self::FEMM,
    assembler::A,
    geom::NodalField{FT},
    afield::NodalField{T},
) where {FEMM<:AbstractFEMM, A<:AbstractSysmatAssembler, FT, T}

Compute the inner-product (Gram) matrix.

inspectintegpoints(self::FEMM,
    geom::NodalField{GFT},
    u::NodalField{FT},
    dT::NodalField{FT},
    felist::AbstractVector{IT},
    inspector::F,
    idat,
    quantity = :Cauchy;
    context...,) where {FEMM<:AbstractFEMM, GFT, IT, FT, F <: Function}

Inspect integration points.

integratefieldfunction(
    self::AbstractFEMM,
    geom::NodalField{GFT},
    afield::FL,
    fh::F;
    initial::R = zero(FT),
    m = -1,
) where {GFT, T, FL<:ElementalField{T}, F<:Function,R}

Integrate a elemental-field function over the discrete manifold.

• afield = ELEMENTAL field to supply the value within the element (one value per element),
• fh = function taking position and an array of field values for the element as arguments, returning value of type R. The function fh must take two arguments, x which is the location, and val which is the value of the field at that location. The rectangular array of field values val has one row, and as many columns as there are degrees of freedom per node.
• m = dimension of the manifold over which to integrate; m < 0 means that the dimension is controlled by the manifold dimension of the elements.

Integrates a function returning a scalar value of type R, which is initialized by initial.

integratefieldfunction(
    self::AbstractFEMM,
    geom::NodalField{GFT},
    afield::FL,
    fh::F;
    initial::R,
    m = -1,
) where {GFT, T, FL<:NodalField{T}, F<:Function,R}

Integrate a nodal-field function over the discrete manifold.

• afield = NODAL field to supply the values at the nodes, which are interpolated to the quadrature points,
• fh = function taking position and an array of field values for the element as arguments, returning value of type R. The function fh must take two arguments, x which is the location, and val which is the value of the field at that location. The rectangular array of field values val has one row, and as many columns as there are degrees of freedom per node.
• m = dimension of the manifold over which to integrate; m < 0 means that the dimension is controlled by the manifold dimension of the elements.

Integrates a function returning a scalar value of type R, which is initialized by initial.

integratefunction(
    self::AbstractFEMM,
    geom::NodalField{GFT},
    fh::F;
    initial::R = zero(typeof(fh(zeros(ndofs(geom), 1)))),
    m = -1,
) where {GFT<:Number, F<:Function, R<:Number}

Integrate a function over the discrete manifold.

Integrate some scalar function over the geometric cells. The function takes a single argument, the position vector.

When the scalar function returns just +1 (such as (x) -> 1.0), the result measures the volume (number of points, length, area, 3-D volume, according to the manifold dimension). When the function returns the mass density, the method measures the mass, when the function returns the x-coordinate equal measure the static moment with respect to the y- axis, and so on.

Example:

Compute the volume of the mesh and then its center of gravity:

V = integratefunction(femm, geom, (x) ->  1.0, 0.0)
Sx = integratefunction(femm, geom, (x) ->  x[1], 0.0)
Sy = integratefunction(femm, geom, (x) ->  x[2], 0.0)
Sz = integratefunction(femm, geom, (x) ->  x[3], 0.0)
CG = vec([Sx Sy Sz]/V)

Compute a moment of inertia of the mesh relative to the origin:

Ixx = integratefunction(femm, geom, (x) ->  x[2]^2 + x[3]^2)

iselementbased(self::FEMM) where {FEMM <: AbstractFEMM}

Is the FEMM element-based? (This will only be false for nodal-integration formulations.)

linform_dot(
    self::FEMM,
    assembler::A,
    geom::NodalField{FT},
    P::NodalField{T},
    f::DC,
    m,
) where {FEMM<:AbstractFEMM, A<:AbstractSysvecAssembler, FT<:Number, T, DC<:DataCache}

Compute the discrete vector implied by the linear form "dot".

\[\int_{V} \mathbf{w} \cdot \mathbf{f} \; \mathrm{d} V\]

Here $\mathbf{w}$ is the test function, $\mathbf{f}$ is a given function (data). Both are assumed to be vectors, even if they are of length 1, representing scalars. The data $\mathbf{f}$ is represented with DataCache.

Arguments

• self = finite element machine;
• assembler = assembler of the global vector;
• geom = geometry field;
• P = nodal field to define the degree of freedom numbers;
• f= data cache, which is called to evaluate the integrand based on the location, the Jacobian matrix, the finite element identifier, and the quadrature point;
• m= manifold dimension, 0= vertex (point), 1= curve, 2= surface, 3= volume. For body loads m is set to 3, for tractions on the surface it is set to 2, and so on.

transferfield!(
    ff::F,
    fensf::FENodeSet{FT},
    fesf::AbstractFESet,
    fc::F,
    fensc::FENodeSet{FT},
    fesc::AbstractFESet,
    geometricaltolerance::FT;
    parametrictolerance::FT = 0.01,
) where {FT<:Number, F<:ElementalField{T}, T}

Transfer an elemental field from a coarse mesh to a finer one.

Arguments

• ff = the fine-mesh field (modified and also returned)
• fensf = finite element node set for the fine-mesh
• fc = the coarse-mesh field
• fensc = finite element node set for the fine-mesh,
• fesc = finite element set for the coarse mesh
• tolerance = tolerance in physical space for searches of the adjacent nodes

Output

Elemental field ff transferred to the fine mesh is output.

transferfield!(
    ff::F,
    fensf::FENodeSet{FT},
    fesf::AbstractFESet,
    fc::F,
    fensc::FENodeSet{FT},
    fesc::AbstractFESet,
    geometricaltolerance::FT;
    parametrictolerance::FT = 0.01,
) where {FT<:Number, F<:NodalField{T}, T}

Transfer a nodal field from a coarse mesh to a finer one.

Arguments

• ff = the fine-mesh field (modified and also returned)
• fensf = finite element node set for the fine-mesh
• fc = the coarse-mesh field
• fensc = finite element node set for the fine-mesh,
• fesc = finite element set for the coarse mesh
• geometricaltolerance = tolerance in physical space for searches of the adjacent nodes
• parametrictolerance = tolerance in parametric space for for check whether node is inside an element

Output

Nodal field ff transferred to the fine mesh is output.

bisect(fun, xl, xu, tolx, tolf)

Implementation of the bisection method.

Tolerance both on x and on f(x) is used.

• fun = function,
• xl= lower value of the bracket,
• xu= upper Value of the bracket,
• tolx= tolerance on the location of the root,
• tolf= tolerance on the function value

bisect(fun, xl, xu, fl, fu, tolx, tolf)

Implementation of the bisection method.

Tolerance both on x and on f(x) is used.

• fun = function,
• xl,xu= lower and upper value of the bracket,
• fl,fu= function value at the lower and upper limit of the bracket.

The true values must have opposite signs (that is they must constitute a bracket). Otherwise this algorithm will fail.

• tolx= tolerance on the location of the root,
• tolf= tolerance on the function value

conjugategradient(A::MT, b::Vector{T}, x0::Vector{T}, maxiter) where {MT, T<:Number}

Compute one or more iterations of the conjugate gradient process.

evalconvergencestudy(modeldatasequence)

Evaluate a convergence study from a model-data sequence.

• modeldatasequence = array of modeldata dictionaries. At least two must be included.

Refer to methods fieldnorm and fielddiffnorm for details on the required keys in the dictionaries.

Output

• elementsizes = element size array,
• errornorms = norms of the error,
• convergencerate = rate of convergence

fielddiffnorm(modeldatacoarse, modeldatafine)

Compute norm of the difference of the fields.

Arguments

• modeldatacoarse, modeldatafine = data dictionaries.

For both the "coarse"- and "fine"-mesh modeldata the data dictionaries need to contain the mandatory keys:

• "fens" = finite element node set
• "regions" = array of regions
• "targetfields" = array of fields, one for each region
• "geom" = geometry field
• "elementsize" = representative element size,
• "geometricaltolerance" = geometrical tolerance (used in field transfer; refer to the documentation of transferfield!)
• "parametrictolerance" = parametric tolerance (used in field transfer; refer to the documentation of transferfield!)

Output

• Norm of the field as floating-point scalar.

fieldnorm(modeldata)

Compute norm of the target field.

Argument

• modeldata = data dictionary, mandatory keys:
  • fens = finite element node set
  • regions = array of regions
  • targetfields = array of fields, one for each region
  • geom = geometry field
  • elementsize = representative element size,

Output

• Norm of the field as floating-point scalar.

matrix_blocked(A, row_nfreedofs, col_nfreedofs = row_nfreedofs)

Partition matrix into blocks.

The function returns the sparse matrix as a named tuple of its constituent blocks. The matrix is assumed to be composed of four blocks

A = [A_ff A_fd
     A_df A_dd]

The named tuple is the value (ff = A_ff, fd = A_fd, df = A_df, dd = A_dd). Index into this named tuple to retrieve the parts of the matrix that you want.

Here f stands for free, and d stands for data (i.e. fixed, prescribed, ...). The size of the ff block is row_nfreedofs, col_nfreedofs. Neither one of the blocks is square, unless row_nfreedofs == col_nfreedofs.

When row_nfreedofs == col_nfreedofs, only the number of rows needs to be given.

Example

Both

K_ff, K_fd = matrix_blocked(K, nfreedofs, nfreedofs)[(:ff, :fd)]
K_ff, K_fd = matrix_blocked(K, nfreedofs)[(:ff, :fd)]

define a square K_ff matrix and, in general a rectangular, matrix K_fd.

This retrieves all four partitions of the matrix

A_ff, A_fd, A_df, A_dd = matrix_blocked(A, nfreedofs)[(:ff, :fd, :df, :dd)]

This retrieves the complete named tuple, and then the matrices can be referenced with a dot syntax.

A_b = matrix_blocked(A, nfreedofs, nfreedofs)
@test size(A_b.ff) == (nfreedofs, nfreedofs)
@test size(A_b.fd) == (nfreedofs, size(A, 1) - nfreedofs)

penaltyebc!(K, F, dofnums, prescribedvalues, penfact)

Apply penalty essential boundary conditions.

Arguments

• K = stiffness matrix
• F = global load vector
• dofnums, prescribedvalues = arrays computed by prescribeddofs()
• penfact = penalty multiplier, in relative terms: how many times the maximum absolute value of the diagonal elements should the penalty term be?

Output

• Updated matrix K and vector F.

qcovariance(ps::VecOrMat{T}, xs::VecOrMat{T}, ys::VecOrMat{T}; ws = nothing) where {T<:Number}

Compute the covariance for two 'functions' given by the arrays xs and ys at the values of the parameter ps. ws is the optional weights vector; if it is not supplied, uniformly distributed default weights are assumed.

Notes:

– The mean is subtracted from both functions. – This function is not particularly efficient: it computes the mean of both functions and it allocates arrays instead of overwriting the contents of the arguments.

qtrap(ps::VecOrMat{T}, xs::VecOrMat{T}) where {T<:Number}

Compute the area under the curve given by a set of parameters along an interval and the values of the 'function' at those parameter values. The parameter values need not be uniformly distributed.

Trapezoidal rule is used to evaluate the integral. The 'function' is assumed to vary linearly inbetween the given points.

qvariance(ps, xs; ws = nothing)

Compute the variance of a function given by the array xs at the values of the parameter ps. ws is the optional weights vector with unit default weights.

richextrapol(solns::T, params::T; lower_conv_rate = 0.001, upper_conv_rate = 10.0) where {T<:AbstractArray{Tn} where {Tn}}

Richardson extrapolation.

Arguments

• solns = array of three solution values
• params = array of values of three parameters for which the solns have been obtained.

The assumption is that the error of the solution is expanded in a Taylor series, and only the first term in the Taylor series is kept. qex - qapprox ~ C param^beta Here qex is the true solution, qapprox is an approximate solution, param is the element size, or the relative element size, in other words the parameter of the extrapolation, and beta is the convergence rate. The constant C is the third unknown quantity in this expansion. If we obtain three successive approximations, we can solve for the three unknown quantities, qex, beta, and C.

It is assumed that the first solution is obtained for the largest value of the extrapolation parameter, while the last solution in the list is obtained for the smallest value of the extrapolation parameter: params[1] > params[2] > params[3]

Output

• solnestim= estimate of the asymptotic solution from the data points in the solns array
• beta= convergence rate
• c = constant in the estimate error=c*h^beta
• maxresidual = maximum residual after equations from which the above quantities were solved (this is a measure of how accurately was the system solved).

richextrapoluniform(solns::T, params::T) where {T<:AbstractArray{Tn} where {Tn}}

Richardson extrapolation.

Argument

• solns = array of solution values
• params = array of values of parameters for which the solns have been obtained. This function is applicable only to fixed (uniform) ratio between the mesh sizes, params[1]/params[2) = params[2)/params[3).

Output

• solnestim= estimate of the asymptotic solution from the data points in the solns array
• beta= convergence rate
• c = constant in the estimate error=c*h^beta
• residual = residual after equations from which the above quantities were solved (this is a measure of how accurately was the system solved).

solve_blocked!(u::AF, K::M, F::V) where {AF<:AbstractField, M<:AbstractMatrix, V<:AbstractVector}

Solve a system of linear algebraic equations.

solve_blocked(A::M, b::VB, x::VX, nfreedofs::IT) where {M<:AbstractMatrix, VB<:AbstractVector, VX<:AbstractVector, IT<:Integer}

Solve a blocked system of linear algebraic equations.

The matrix is blocked as

A = [A_ff A_fd
     A_df A_dd]

and the solution and the righthand side vector are blocked accordingly

x = [x_f
     x_d]

and

b = [b_f
     b_d]

Above, b_f andxdare known,xf(solution) andb_d` (reactions) need to be computed.

vector_blocked(V, row_nfreedofs, which = (:all, ))

Partition vector into two pieces.

The vector is composed of two blocks

V = [V_f
     V_d]

which are returned as a named tuple (f = V_f, d = V_d).

massdensity(self::AbstractMat)

Return mass density.