The adaptivity module provides a framework to work with adapted (refined/coarsened/mixed) meshes.
It provides
A generic interface to represent adapted meshes and a set of tools to work with Finite Element spaces defined on them. In particular, moving CellFields between parent and child meshes.
Particular implementations for conformally refining/coarsening 2D/3D meshes using several well-known strategies. In particular, Red-Green refinement and longest-edge bisection.
Glue containing the map between two nested triangulations. The contained datastructures will depend on the type of glue. There are two types of AdaptivityGlue:
RefinementGlue :: All cells in the new mesh are children of cells in the old mesh. I.e given a new cell, it is possible to find a single old cell containing it (the new cell might be exactly the old cell if no refinement).
MixedGlue :: Some cells in the new mesh are children of cells in the old mesh, while others are parents of cells in the old mesh.
Contains:
n2o_faces_map :: Given a new face gid, returns
if fine, the gid of the old face containing it.
if coarse, the gids of its children (in child order)
n2o_cell_to_child_id :: Given a new cell gid, returns
if fine, the local child id within the (old) coarse cell containing it.
if coarse, a list of local child ids of the (old) cells containing it.
refinement_rules :: Array conatining the RefinementRule used for each coarse cell.
DiscreteModel created by refining/coarsening another DiscreteModel.
The refinement/coarsening hierarchy can be traced backwards by following the parent pointer chain. This allows the transfer of dofs between FESpaces defined on this model and its ancestors.
The module provides a refine method for UnstructuredDiscreteModel. The method takes a string refinement_method that determines the refinement strategy to be used. The following strategies are available:
"red_green" :: Red-Green refinement, default.
"nvb" :: Longest-edge bisection (only for meshes of TRIangles)
"barycentric" :: Barycentric refinement (only for meshes of TRIangles)
"simplexify" :: Simplexify refinement. Same resulting mesh as the simplexify method, but keeps track of the parent-child relationships.
Additionally, the method takes a kwarg cells_to_refine that determines which cells will be refined. Possible input types are:
Nothing :: All cells get refined.
AbstractArray{<:Bool} of size num_cells(model) :: Only cells such that cells_to_refine[iC] == true get refined.
AbstractArray{<:Integer} :: Cells for which gid ∈ cells_to_refine get refined
The algorithms try to respect the cells_to_refine input as much as possible, but some additional cells might get refined in order to guarantee that the mesh remains conforming.
function refine(model::UnstructuredDiscreteModel;refinement_method="red_green",kwargs...)
+Gridap.Adaptivity · Gridap.jl
The adaptivity module provides a framework to work with adapted (refined/coarsened/mixed) meshes.
It provides
A generic interface to represent adapted meshes and a set of tools to work with Finite Element spaces defined on them. In particular, moving CellFields between parent and child meshes.
Particular implementations for conformally refining/coarsening 2D/3D meshes using several well-known strategies. In particular, Red-Green refinement and longest-edge bisection.
Glue containing the map between two nested triangulations. The contained datastructures will depend on the type of glue. There are two types of AdaptivityGlue:
RefinementGlue :: All cells in the new mesh are children of cells in the old mesh. I.e given a new cell, it is possible to find a single old cell containing it (the new cell might be exactly the old cell if no refinement).
MixedGlue :: Some cells in the new mesh are children of cells in the old mesh, while others are parents of cells in the old mesh.
Contains:
n2o_faces_map :: Given a new face gid, returns
if fine, the gid of the old face containing it.
if coarse, the gids of its children (in child order)
n2o_cell_to_child_id :: Given a new cell gid, returns
if fine, the local child id within the (old) coarse cell containing it.
if coarse, a list of local child ids of the (old) cells containing it.
refinement_rules :: Array conatining the RefinementRule used for each coarse cell.
DiscreteModel created by refining/coarsening another DiscreteModel.
The refinement/coarsening hierarchy can be traced backwards by following the parent pointer chain. This allows the transfer of dofs between FESpaces defined on this model and its ancestors.
The module provides a refine method for UnstructuredDiscreteModel. The method takes a string refinement_method that determines the refinement strategy to be used. The following strategies are available:
"red_green" :: Red-Green refinement, default.
"nvb" :: Longest-edge bisection (only for meshes of TRIangles)
"barycentric" :: Barycentric refinement (only for meshes of TRIangles)
"simplexify" :: Simplexify refinement. Same resulting mesh as the simplexify method, but keeps track of the parent-child relationships.
Additionally, the method takes a kwarg cells_to_refine that determines which cells will be refined. Possible input types are:
Nothing :: All cells get refined.
AbstractArray{<:Bool} of size num_cells(model) :: Only cells such that cells_to_refine[iC] == true get refined.
AbstractArray{<:Integer} :: Cells for which gid ∈ cells_to_refine get refined
The algorithms try to respect the cells_to_refine input as much as possible, but some additional cells might get refined in order to guarantee that the mesh remains conforming.
function refine(model::UnstructuredDiscreteModel;refinement_method="red_green",kwargs...)
[...]
end
The module provides a refine method for CartesianDiscreteModel. The method takes a Tuple of size Dc (the dimension of the model cells) that will determine how many times cells will be refined in each direction. For example, for a 2D model, refine(model,(2,3)) will refine each QUAD cell into a 2x3 grid of cells.
function refine(model::CartesianDiscreteModel{Dc}, cell_partition::Tuple) where Dc
[...]
- end
The module also provides support for macro finite-elements. From an abstract point of view, a macro finite-element is a finite-element defined on a refined polytope, where polynomial basis are defined on each of the subcells (creating a broken piece-wise polynomial space on the original polytope). From Gridap's point of view, a macro finite-element is a ReferenceFE defined on a RefinementRule from an array of ReferenceFEs defined on the subcells.
Although there are countless combinations, here are two possible applications:
Linearized High-Order Lagrangian FESpaces: These are spaces which have the same DoFs as a high-order Lagrangian space, but where the basis functions are linear on each subcell.
Barycentrically-refined elements for Stokes-like problems: These are spaces where the basis functions for velocity are defined on the barycentrically-refined mesh, whereas the basis functions for pressure are defined on the original cells. This allows for exact so-called Stokes sequences (see here).
Most of the tools provided by this module are showcased in the tests of the module itself, as well as the following tutorial (coming soon).
However, we want to stress a couple of key performance-critical points:
The refining/coarsening routines are not optimized for performance. In particular, they are not parallelized. If you require an optimized/parallel implementation, please consider leveraging specialised meshing libraries. For instance, we provide an implementation of refine/coarsen using P4est in the GridapP4est.jl library.
Although the toolbox allows you to evaluate CellFields defined on both fine/coarse meshes on their parent/children mesh, both directions of evaluation are not equivalent. As a user, you should always try to evaluate/integrate on the finest mesh for maximal performance. Evaluating a fine CellField on a coarse mesh relies on local tree searches, and is therefore a very expensive operation that should be avoided whenever possible.
Given a RefinementRule, the library provides a set of methods to compute the mappings between parent (coarse) face ids and child (fine) face ids (and vice-versa).
The most basic information (that can directly be hardcoded in the RefinementRule for performance) are the mappings between parent face ids and child face ids. These are provided by:
Given a RefinementRule, returns for each parent/coarse face the child/fine faces of the same dimension that it contains. Therefore, only fine faces at the coarse cell boundary are listed in the returned structure.
Returns: [Face dimension][Coarse Face id] -> [Fine faces]
On top of these two basic mappings, a whole plethora of additional topological mappings can be computed. These first set of routines extend the ReferenceFEs API to provide information on the face-to-node mappings and permutations:
Given a polytope p, returns a vector of vectors containing all admissible permutations of the polytope vertices. An admissible permutation is one such that, if the vertices of the polytope are re-labeled according to this permutation, the resulting polytope preserves the shape of the original one.
The module also provides support for macro finite-elements. From an abstract point of view, a macro finite-element is a finite-element defined on a refined polytope, where polynomial basis are defined on each of the subcells (creating a broken piece-wise polynomial space on the original polytope). From Gridap's point of view, a macro finite-element is a ReferenceFE defined on a RefinementRule from an array of ReferenceFEs defined on the subcells.
Although there are countless combinations, here are two possible applications:
Linearized High-Order Lagrangian FESpaces: These are spaces which have the same DoFs as a high-order Lagrangian space, but where the basis functions are linear on each subcell.
Barycentrically-refined elements for Stokes-like problems: These are spaces where the basis functions for velocity are defined on the barycentrically-refined mesh, whereas the basis functions for pressure are defined on the original cells. This allows for exact so-called Stokes sequences (see here).
Most of the tools provided by this module are showcased in the tests of the module itself, as well as the following tutorial (coming soon).
However, we want to stress a couple of key performance-critical points:
The refining/coarsening routines are not optimized for performance. In particular, they are not parallelized. If you require an optimized/parallel implementation, please consider leveraging specialised meshing libraries. For instance, we provide an implementation of refine/coarsen using P4est in the GridapP4est.jl library.
Although the toolbox allows you to evaluate CellFields defined on both fine/coarse meshes on their parent/children mesh, both directions of evaluation are not equivalent. As a user, you should always try to evaluate/integrate on the finest mesh for maximal performance. Evaluating a fine CellField on a coarse mesh relies on local tree searches, and is therefore a very expensive operation that should be avoided whenever possible.
Given a RefinementRule, the library provides a set of methods to compute the mappings between parent (coarse) face ids and child (fine) face ids (and vice-versa).
The most basic information (that can directly be hardcoded in the RefinementRule for performance) are the mappings between parent face ids and child face ids. These are provided by:
Given a RefinementRule, returns for each parent/coarse face the child/fine faces of the same dimension that it contains. Therefore, only fine faces at the coarse cell boundary are listed in the returned structure.
Returns: [Face dimension][Coarse Face id] -> [Fine faces]
On top of these two basic mappings, a whole plethora of additional topological mappings can be computed. These first set of routines extend the ReferenceFEs API to provide information on the face-to-node mappings and permutations:
Given a polytope p, returns a vector of vectors containing all admissible permutations of the polytope vertices. An admissible permutation is one such that, if the vertices of the polytope are re-labeled according to this permutation, the resulting polytope preserves the shape of the original one.
The first admissible permutation for a segment is [1,2],i.e., the identity. The second one is [2,1], i.e., the first vertex is relabeled as 2 and the second vertex is relabeled as 1.
The first admissible permutation for a segment is [1,2],i.e., the identity. The second one is [2,1], i.e., the first vertex is relabeled as 2 and the second vertex is relabeled as 1.
aggregate_cface_to_own_fface_data(
rr::RefinementRule,
cface_to_own_fface_to_data :: AbstractVector{<:AbstractVector{T}}
-) where T
Given a RefinementRule, and a data structure cface_to_own_fface_to_data that contains data for each child/fine face owned by each parent/coarse face, returns a data structure cface_to_fface_to_data that contains the data for each child/fine face contained in the closure of each parent/coarse face (i.e the fine faces are owned and not owned).
The implementation makes sure that the resulting data structure is ordered according to the fine face numbering in get_cface_to_ffaces(rrule) (which in turn is by increasing fine face id).
Given a RefinementRule of dimension Dc and a Dc-Tuple fine_orders of approximation orders, returns a map between the fine nodal dofs of order fine_orders in the reference grid and the coarse nodal dofs of order 2⋅fine_orders in the coarse parent cell.
The result is given for each coarse/parent face of dimension D as a list of the corresponding fine dof lids, i.e
When a cell is refined, we need to be able to evaluate the fields defined on the children cells on the parent cell. To do so, we bundle the fields defined on the children cells into a new type of Field called FineToCoarseField. When evaluated on a Point, a FineToCoarseField will select the child cell that contains the Point and evaluate the mapped point on the corresponding child field.
Given a RefinementRule, and a data structure cface_to_own_fface_to_data that contains data for each child/fine face owned by each parent/coarse face, returns a data structure cface_to_fface_to_data that contains the data for each child/fine face contained in the closure of each parent/coarse face (i.e the fine faces are owned and not owned).
The implementation makes sure that the resulting data structure is ordered according to the fine face numbering in get_cface_to_ffaces(rrule) (which in turn is by increasing fine face id).
Given a RefinementRule of dimension Dc and a Dc-Tuple fine_orders of approximation orders, returns a map between the fine nodal dofs of order fine_orders in the reference grid and the coarse nodal dofs of order 2⋅fine_orders in the coarse parent cell.
The result is given for each coarse/parent face of dimension D as a list of the corresponding fine dof lids, i.e
When a cell is refined, we need to be able to evaluate the fields defined on the children cells on the parent cell. To do so, we bundle the fields defined on the children cells into a new type of Field called FineToCoarseField. When evaluated on a Point, a FineToCoarseField will select the child cell that contains the Point and evaluate the mapped point on the corresponding child field.
Given a domain and a non-overlapping refined cover, a FineToCoarseField is a Field defined in the domain and constructed by a set of fields defined on the subparts of the covering partition. The refined cover is represented by a RefinementRule.
Parameters:
rrule: Refinement rule representing the covering partition.
fine_fields: Fields defined on the subcells of the covering partition. To accomodate the case where not all subcells are present (e.g in distributed), we allow for length(fine_fields) != num_subcells(rrule).
id_map: Mapping from the subcell ids to the indices of the fine_fields array. If length(fine_fields) == num_subcells(rrule), this is the identity map. Otherwise, the id_map is used to map the subcell ids to the indices of the fine_fields array.
This document was generated with Documenter.jl version 1.8.0 on Thursday 28 November 2024. Using Julia version 1.10.7.
+end
Given a domain and a non-overlapping refined cover, a FineToCoarseField is a Field defined in the domain and constructed by a set of fields defined on the subparts of the covering partition. The refined cover is represented by a RefinementRule.
Parameters:
rrule: Refinement rule representing the covering partition.
fine_fields: Fields defined on the subcells of the covering partition. To accomodate the case where not all subcells are present (e.g in distributed), we allow for length(fine_fields) != num_subcells(rrule).
id_map: Mapping from the subcell ids to the indices of the fine_fields array. If length(fine_fields) == num_subcells(rrule), this is the identity map. Otherwise, the id_map is used to map the subcell ids to the indices of the fine_fields array.
The cache generated when using this solver has a field result that hosts the result object generated by the underlying nlsolve function. It corresponds to the most latest solve.
The cache generated when using this solver has a field result that hosts the result object generated by the underlying nlsolve function. It corresponds to the most latest solve.
Type representing an array with a reduced set of values. The array is represented by a short array of values, namely the field values, and a large array of indices, namely the field ptrs. The i-th component of the resulting array is defined as values[ptrs[i]]. The type parameters A, and P are restricted to be array types by the inner constructor of this struct.
Subtype of AbstractArray which is the result of lazy_map. It represents the result of lazy_mapping a Map to a set of arrays that contain the mapping arguments. This struct makes use of the cache provided by the mapping in order to compute its indices (thus allowing to prevent allocation). The array is lazy, i.e., the values are only computed on demand. It extends the AbstractArray API with two methods:
Abstract type representing a function (mapping) that provides a cache and an in-place evaluation for performance. This is the type to be used in the lazy_map function.
Derived types must implement the following method:
The mapping interface can be tested with the test_map function.
Note that most of the functionality implemented in terms of this interface relies in duck typing. That is, it is not strictly needed to work with types that inherit from Map. This is specially useful in order to accommodate existing types into this framework without the need to implement a wrapper type that inherits from Map. For instance, a default implementation is available for Function objects. However, we recommend that new types inherit from Map.
Returns the map that results after applying an operation f over a set of map(s) args. That is Operation(f)(args)(x...) is formally defined as f(map(a->a(x...),args)...).
Example
using Gridap.Arrays
+end
Type representing an array with a reduced set of values. The array is represented by a short array of values, namely the field values, and a large array of indices, namely the field ptrs. The i-th component of the resulting array is defined as values[ptrs[i]]. The type parameters A, and P are restricted to be array types by the inner constructor of this struct.
Subtype of AbstractArray which is the result of lazy_map. It represents the result of lazy_mapping a Map to a set of arrays that contain the mapping arguments. This struct makes use of the cache provided by the mapping in order to compute its indices (thus allowing to prevent allocation). The array is lazy, i.e., the values are only computed on demand. It extends the AbstractArray API with two methods:
Abstract type representing a function (mapping) that provides a cache and an in-place evaluation for performance. This is the type to be used in the lazy_map function.
Derived types must implement the following method:
The mapping interface can be tested with the test_map function.
Note that most of the functionality implemented in terms of this interface relies in duck typing. That is, it is not strictly needed to work with types that inherit from Map. This is specially useful in order to accommodate existing types into this framework without the need to implement a wrapper type that inherits from Map. For instance, a default implementation is available for Function objects. However, we recommend that new types inherit from Map.
Returns the map that results after applying an operation f over a set of map(s) args. That is Operation(f)(args)(x...) is formally defined as f(map(a->a(x...),args)...).
Returns a mapping that represents the result of applying the function f to the arguments in the tuple args. That is, OperationMap(f,args)(x...) is formally defined as f(map(a->a(x...),args)...)
Returns a mapping that represents the result of applying the function f to the arguments in the tuple args. That is, OperationMap(f,args)(x...) is formally defined as f(map(a->a(x...),args)...)
Returns a cache object to be used in the getindex! function. It defaults to
array_cache(a::T) where T = nothing
for types T such that uses_hash(T) == Val(false), and
function array_cache(a::T) where T
hash = Dict{UInt,Any}()
array_cache(hash,a)
end
for types T such that uses_hash(T) == Val(true), see the uses_hash function. In the later case, the type T should implement the following signature:
array_cache(hash::Dict,a::AbstractArray)
where we pass a dictionary (i.e., a hash table) in the first argument. This hash table can be used to test if the object a has already built a cache and re-use it as follows
id = objectid(a)
@@ -54,10 +54,10 @@
else
cache = ... # Build a new cache depending on your needs
hash[id] = cache # Register the cache in the hash table
-end
This mechanism is needed, e.g., to re-use intermediate results in complex lazy operation trees. In multi-threading computations, a different hash table per thread has to be used in order to avoid race conditions.
Applies the mapping f at the arguments x... using the scratch data provided in the given cache object. The cache object is built with the return_cache function using arguments of the same type as in x. In general, the returned value y can share some part of its state with the cache object. If the result of two or more calls to this function need to be accessed simultaneously (e.g., in multi-threading), create and use several cache objects (e.g., one cache per thread).
Returns the item of the array a associated with index i by (possibly) using the scratch data passed in the cache object.
It defaults to
getindex!(cache,a::AbstractArray,i...) = a[i...]
As for standard Julia arrays, the user needs to implement only one of the following signatures depending on the IndexStyle of the array.
getindex!(cache,a::AbstractArray,i::Integer)
+end
This mechanism is needed, e.g., to re-use intermediate results in complex lazy operation trees. In multi-threading computations, a different hash table per thread has to be used in order to avoid race conditions.
Applies the mapping f at the arguments x... using the scratch data provided in the given cache object. The cache object is built with the return_cache function using arguments of the same type as in x. In general, the returned value y can share some part of its state with the cache object. If the result of two or more calls to this function need to be accessed simultaneously (e.g., in multi-threading), create and use several cache objects (e.g., one cache per thread).
Applies the Map (or Function) f to the entries of the arrays in a (see the definition of Map).
The resulting array r is such that r[i] equals to evaluate(f,ai...) where ai is the tuple containing the i-th entry of the arrays in a (see function evaluate for more details). In other words, the resulting array is numerically equivalent to:
Applies the Map (or Function) f to the entries of the arrays in a (see the definition of Map).
The resulting array r is such that r[i] equals to evaluate(f,ai...) where ai is the tuple containing the i-th entry of the arrays in a (see function evaluate for more details). In other words, the resulting array is numerically equivalent to:
map( (x...)->evaluate(f,x...), a...)
Examples
Using a function as mapping
using Gridap.Arrays
a = collect(0:5)
b = collect(10:15)
@@ -101,9 +101,9 @@
println(c)
# output
-[10, 12, 14, 16, 18, 20]
Returns the cache needed to lazy_map mapping f with arguments of the same type as the objects in x. This function returns nothing by default, i.e., no cache.
Function used to test if the mapping f has been implemented correctly. f is a Map sub-type, x is a tuple in the domain of the mapping and y is the expected result. Function cmp is used to compare the computed result with the expected one. The checks are done with the @test macro.
The default implementation of this function is testargs(f,x...) = x. One can overload it in order to use lazy_map with 0-length array and maps with non-trivial domains.
Returns the cache needed to lazy_map mapping f with arguments of the same type as the objects in x. This function returns nothing by default, i.e., no cache.
Function used to test if the mapping f has been implemented correctly. f is a Map sub-type, x is a tuple in the domain of the mapping and y is the expected result. Function cmp is used to compare the computed result with the expected one. The checks are done with the @test macro.
The default implementation of this function is testargs(f,x...) = x. One can overload it in order to use lazy_map with 0-length array and maps with non-trivial domains.
Returns an arbitrary instance of eltype(a). The default returned value is the first entry in the array if length(a)>0 and testvalue(eltype(a)) if length(a)==0 See the testvalue function.
Examples
using Gridap.Arrays
a = collect(3:10)
@@ -116,4 +116,4 @@
# output
(3, 0)
-
Returns an arbitrary instance of type T. It defaults to zero(T) for non-array types and to an empty array for array types. It can be overloaded for new types T if zero(T) does not makes sense. This function is used to compute testitem for 0-length arrays.
This function is used to specify if the type T uses the hash-based mechanism to reuse caches. It should return either Val(true) or Val(false). It defaults to
uses_hash(::Type{<:AbstractArray}) = Val(false)
Once this function is defined for the type T it can also be called on instances of T.
Returns an arbitrary instance of type T. It defaults to zero(T) for non-array types and to an empty array for array types. It can be overloaded for new types T if zero(T) does not makes sense. This function is used to compute testitem for 0-length arrays.
This function is used to specify if the type T uses the hash-based mechanism to reuse caches. It should return either Val(true) or Val(false). It defaults to
uses_hash(::Type{<:AbstractArray}) = Val(false)
Once this function is defined for the type T it can also be called on instances of T.
Data associated with the cells of a Triangulation. CellDatum objects behave as if they are defined in the physical space of the triangulation. But in some cases they are implemented as reference quantities plus some transformation to the physical domain.
SkeletonCellFieldPair is a special construct for allowing uh.plus and uh.minus to be two different CellFields. In particular, it is useful when we need one of the CellFields to be the dualized version of the other for performing ForwardDiff AD of a Skeleton integration DomainContribution wrt to the degrees of freedom of the CellField, plus and minus sensitivities done separately, so as to restrict the interaction between the dual numbers.
It takes in two CellFields and stores plus version of CellFieldAt of the first CellField and minus version of CellFieldAt of the second the CellField. SkeletonCellFieldPair is associated with same triangulation as that of the CellFields (we check if the triangulations of both CellFields match)
SkeletonCellFieldPair is an internal convenience artifact/construct to aid in dualizing plus and minus side around a Skeleton face separately to perform the sensitivity of degrees of freedom of cells sharing the Skeleton face, without the interaction dual numbers of the two cells. The user doesn't have to deal with this construct anywhere when performing AD of functionals involving integration over Skeleton faces using the public API.
Calculate distance from point x to the polytope. The polytope is given by its type and by the inverse cell map, i.e. by the map from the physical to the reference space.
Positive distances are outside the polytope, negative distances are inside the polytope.
The distance is measured in an unspecified norm, currently the L∞ norm.
Data associated with the cells of a Triangulation. CellDatum objects behave as if they are defined in the physical space of the triangulation. But in some cases they are implemented as reference quantities plus some transformation to the physical domain.
SkeletonCellFieldPair is a special construct for allowing uh.plus and uh.minus to be two different CellFields. In particular, it is useful when we need one of the CellFields to be the dualized version of the other for performing ForwardDiff AD of a Skeleton integration DomainContribution wrt to the degrees of freedom of the CellField, plus and minus sensitivities done separately, so as to restrict the interaction between the dual numbers.
It takes in two CellFields and stores plus version of CellFieldAt of the first CellField and minus version of CellFieldAt of the second the CellField. SkeletonCellFieldPair is associated with same triangulation as that of the CellFields (we check if the triangulations of both CellFields match)
SkeletonCellFieldPair is an internal convenience artifact/construct to aid in dualizing plus and minus side around a Skeleton face separately to perform the sensitivity of degrees of freedom of cells sharing the Skeleton face, without the interaction dual numbers of the two cells. The user doesn't have to deal with this construct anywhere when performing AD of functionals involving integration over Skeleton faces using the public API.
Calculate distance from point x to the polytope. The polytope is given by its type and by the inverse cell map, i.e. by the map from the physical to the reference space.
Positive distances are outside the polytope, negative distances are inside the polytope.
The distance is measured in an unspecified norm, currently the L∞ norm.
Minimum data required to describe dof ownership. At this moment, the cell-wise ownership is compressed on cell types. This can be relaxed in the future, to have an arbitrary cell-wise dof ownership.
Minimum data required to build a conforming FE space. At this moment, the some cell-wise info is compressed on cell types. This can be relaxed in the future, and have an arbitrary cell-wise data.
Minimum data required to describe dof ownership. At this moment, the cell-wise ownership is compressed on cell types. This can be relaxed in the future, to have an arbitrary cell-wise dof ownership.
Minimum data required to build a conforming FE space. At this moment, the some cell-wise info is compressed on cell types. This can be relaxed in the future, and have an arbitrary cell-wise data.
Given a Discrete Model and a reffe, builds a new grid in which the geometrical map is a FEFunction. This is useful when considering geometrical maps that are the result of a FE problem (mesh displacement).
Given a Discrete Model and a reffe, builds a new grid in which the geometrical map is a FEFunction. This is useful when considering geometrical maps that are the result of a FE problem (mesh displacement).
Compute the residual and jacobian of an operator op at a given point u. Depending on the nature of op the point u can either be a plain array or a FEFunction.
This function changes the state of the input and can render it in a corrupted state. It is recommended to rewrite the input uh with the output as illustrated to prevent any issue. If cache===nothing, then it creates a new cache object.
This function changes the state of the input and can render it in a corrupted state. It is recommended to rewrite the input uh with the output as illustrated to prevent any issue.
Compute the residual and jacobian of an operator op at a given point u. Depending on the nature of op the point u can either be a plain array or a FEFunction.
This function changes the state of the input and can render it in a corrupted state. It is recommended to rewrite the input uh with the output as illustrated to prevent any issue. If cache===nothing, then it creates a new cache object.
This function changes the state of the input and can render it in a corrupted state. It is recommended to rewrite the input uh with the output as illustrated to prevent any issue.
Return an "algebraic view" of an operator. Algebraic means, that the resulting operator acts on plain arrays, instead of FEFunctions. This can be useful for solving with external tools like NLsolve.jl. See also FEOperator.
like interpolate, but also compute new degrees of freedom for the dirichlet component. The resulting FEFunction does not necessary belongs to the underlying space
Return an "algebraic view" of an operator. Algebraic means, that the resulting operator acts on plain arrays, instead of FEFunctions. This can be useful for solving with external tools like NLsolve.jl. See also FEOperator.
like interpolate, but also compute new degrees of freedom for the dirichlet component. The resulting FEFunction does not necessary belongs to the underlying space
Abstract type representing a physical (scalar, vector, or tensor) field. The domain is a Point and the range a scalar (i.e., a sub-type of Julia Number), a VectorValue, or a TensorValue.
These different cases are distinguished by the return value obtained when evaluating them. E.g., a physical field returns a vector of values when evaluated at a vector of points, and a basis of nf fields returns a 2d matrix (np x nf) when evaluated at a vector of np points.
The following functions (i.e., the Map API) need to be overloaded:
These four methods are only designed to be called by the default implementation of field_gradient(f) and thus cannot be assumed that they are available for an arbitrary field. For this reason, these functions are not exported. The general way of evaluating a gradient of a field is to build the gradient with gradient(f) and evaluating the resulting object. For evaluating the hessian, use two times gradient.
For performance, the user can also consider a vectorised version of the Field API that evaluates the field in a vector of points (instead of only one point). E.g., the evaluate! function for a vector of points returns a vector of scalar, vector or tensor values.
Given a matrix np x nf1 x nf2 result of the evaluation of a field vector on a vector of points, it returns an array in which the field axes (second and third axes) are permuted. It is equivalent as Base.permutedims(A,(1,3,2)) but more performant, since it does not involve allocations.
Implementation of return_cache for a array of Field.
If the field vector has length nf and it is evaluated in one point, it returns an nf vector with the result. If the same array is applied to a vector of np points, it returns a matrix np x nf.
Abstract type representing a physical (scalar, vector, or tensor) field. The domain is a Point and the range a scalar (i.e., a sub-type of Julia Number), a VectorValue, or a TensorValue.
These different cases are distinguished by the return value obtained when evaluating them. E.g., a physical field returns a vector of values when evaluated at a vector of points, and a basis of nf fields returns a 2d matrix (np x nf) when evaluated at a vector of np points.
The following functions (i.e., the Map API) need to be overloaded:
These four methods are only designed to be called by the default implementation of field_gradient(f) and thus cannot be assumed that they are available for an arbitrary field. For this reason, these functions are not exported. The general way of evaluating a gradient of a field is to build the gradient with gradient(f) and evaluating the resulting object. For evaluating the hessian, use two times gradient.
For performance, the user can also consider a vectorised version of the Field API that evaluates the field in a vector of points (instead of only one point). E.g., the evaluate! function for a vector of points returns a vector of scalar, vector or tensor values.
Given a matrix np x nf1 x nf2 result of the evaluation of a field vector on a vector of points, it returns an array in which the field axes (second and third axes) are permuted. It is equivalent as Base.permutedims(A,(1,3,2)) but more performant, since it does not involve allocations.
Implementation of return_cache for a array of Field.
If the field vector has length nf and it is evaluated in one point, it returns an nf vector with the result. If the same array is applied to a vector of np points, it returns a matrix np x nf.
Function used to test the field interface. v is an array containing the expected result of evaluating the field f at the point or vector of points x. The comparison is performed using the cmp function. For fields objects that support the gradient function, the keyword argument grad can be used. It should contain the result of evaluating gradient(f) at x. Idem for gradgrad. The checks are performed with the @test macro.
This document was generated with Documenter.jl version 1.8.0 on Thursday 28 November 2024. Using Julia version 1.10.7.
+ gradgrad=nothing)
Function used to test the field interface. v is an array containing the expected result of evaluating the field f at the point or vector of points x. The comparison is performed using the cmp function. For fields objects that support the gradient function, the keyword argument grad can be used. It should contain the result of evaluating gradient(f) at x. Idem for gradgrad. The checks are performed with the @test macro.
Builds a CartesianDiscreteModel object which represents a subgrid of a (larger) grid represented by desc. This subgrid is described by its D-dimensional minimum (cmin) and maximum (cmax) CartesianIndex identifiers.
Builds a CartesianDiscreteModel object which represents a subgrid of a (larger) grid represented by desc. This subgrid is described by its D-dimensional minimum (cmin) and maximum (cmax) CartesianIndex identifiers.
Abstract type holding information about a physical grid, the underlying grid topology, and a labeling of the grid faces. This is the information that typically provides a mesh generator, and it is what one needs to perform a simulation.
The DiscreteModel interface is defined by overloading the methods:
Abstract type holding information about a physical grid, the underlying grid topology, and a labeling of the grid faces. This is the information that typically provides a mesh generator, and it is what one needs to perform a simulation.
The DiscreteModel interface is defined by overloading the methods:
function UnstructuredGrid(
node_coordinates::Vector{Point{Dp,Tp}},
cell_node_ids::Table{Ti},
reffes::Vector{<:LagrangianRefFE{Dc}},
cell_types::Vector,
orientation_style::OrientationStyle=NonOriented()) where {Dc,Dp,Tp,Ti}
-end
Given an n-cube type ExtrusionPolytope{D}, returns V=Vector{CartesianIndex} with as many entries as the number of faces in the boundary of the Polytope. For an entry facelid in this vector, V[facelid] returns what has to be added to the CartesianIndex of a cell in order to obtain the CartesianIndex of the cell neighbour of K across the face F with local ID face_lid.
Given an n-cube type ExtrusionPolytope{D}, returns V=Vector{CartesianIndex} with as many entries as the number of faces in the boundary of the Polytope. For an entry facelid in this vector, V[facelid] returns what has to be added to the CartesianIndex of a cell in order to obtain the CartesianIndex of the cell neighbour of K across the face F with local ID face_lid.
If possible, returns a Triangulation to which CellDatum objects can be transferred from trian1 and trian2. Can be trian1, trian2 or a new Triangulation.
If possible, returns a Triangulation to which CellDatum objects can be transferred from trian1 and trian2. Can be trian1, trian2 or a new Triangulation.
The first of the given tags appearing in the face is taken. If there is no tag on a face, this face will have a value equal to UNSET. If not tag or tags are provided, all the tags in the model are considered
The first of the given tags appearing in the face is taken. If there is no tag on a face, this face will have a value equal to UNSET. If not tag or tags are provided, all the tags in the model are considered
Sets the execution mode to either "debug" or "performance", which controls the behavior of the @check macro within the Gridap package.
Debug mode (default): The @check macro will be active, which activates consistency checks within the library. This mode is recommended for development and debugging purposes.
Performance mode: The @check macro will be deactivated. This mode is recommended for production runs, where no errors are expected.
Pre-defined functions set_debug_mode and set_performance_mode are also available. Feature only available in Julia 1.6 and later due to restrictions from Preferences.jl.
Returns a tuple of length D that contains D times the object v. In contrast to tuple(fill(v,D)...) which returns the same result, this function is type-stable.
Sets the execution mode to either "debug" or "performance", which controls the behavior of the @check macro within the Gridap package.
Debug mode (default): The @check macro will be active, which activates consistency checks within the library. This mode is recommended for development and debugging purposes.
Performance mode: The @check macro will be deactivated. This mode is recommended for production runs, where no errors are expected.
Pre-defined functions set_debug_mode and set_performance_mode are also available. Feature only available in Julia 1.6 and later due to restrictions from Preferences.jl.
Returns a tuple of length D that contains D times the object v. In contrast to tuple(fill(v,D)...) which returns the same result, this function is type-stable.
Check validity of a dictionary dict for an object of type T. It runs successfully if the dictionary is valid for a particular type or throws an error in any other case.
De-serialize an object of type T from the dictionary dict. Values stored into this dictionary must be of any native Julia data type (Real, Integer, String, etc.) Dictionary keys are Symbols.
function from_jld2_file(filename::AbstractString,dataset::AbstractString="data")
Loads an object from a JLD2 file given its filename and, optionally, the dataset name. The dataset specifies the name and the root location of the data inside the generated JLD2 file.
function from_jld2_file(::Type{T},filename::AbstractString,dataset::AbstractString="data") where T
Loads an object from a JLD2 file given its filename and, optionally, the dataset name. The dataset specifies the name and the root location of the data inside the generated JLD2 file. Checks if the returned object is of the expected Type{T}, if not return error.
Serialize object into a dictionary of type Dict{Symbol,Any}. Values stored into this dictionary must be of any native Julia data type (Real, Integer, String, etc.) Dictionary keys are Symbols.
function to_jld2_file(object,filename::AbstractString,dataset::AbstractString="data")
Stores an object to a JLD2 file given its filename and, optionally, the dataset name. The dataset specifies the name and the root location of the data inside the JLD2 file.
Check validity of a dictionary dict for an object of type T. It runs successfully if the dictionary is valid for a particular type or throws an error in any other case.
De-serialize an object of type T from the dictionary dict. Values stored into this dictionary must be of any native Julia data type (Real, Integer, String, etc.) Dictionary keys are Symbols.
function from_jld2_file(filename::AbstractString,dataset::AbstractString="data")
Loads an object from a JLD2 file given its filename and, optionally, the dataset name. The dataset specifies the name and the root location of the data inside the generated JLD2 file.
function from_jld2_file(::Type{T},filename::AbstractString,dataset::AbstractString="data") where T
Loads an object from a JLD2 file given its filename and, optionally, the dataset name. The dataset specifies the name and the root location of the data inside the generated JLD2 file. Checks if the returned object is of the expected Type{T}, if not return error.
Serialize object into a dictionary of type Dict{Symbol,Any}. Values stored into this dictionary must be of any native Julia data type (Real, Integer, String, etc.) Dictionary keys are Symbols.
function to_jld2_file(object,filename::AbstractString,dataset::AbstractString="data")
Stores an object to a JLD2 file given its filename and, optionally, the dataset name. The dataset specifies the name and the root location of the data inside the JLD2 file.
Similar to ConsecutiveMultiFieldStyle, but we keep the original DoF ids of the individual spaces for better block assembly (see BlockSparseMatrixAssembler). Takes three parameters:
NB: Number of assembly blocks
SB: Size of each assembly block, as a Tuple.
P : Permutation of the variables of the multifield space when assembling, as a Tuple.
Similar to ConsecutiveMultiFieldStyle, but we keep the original DoF ids of the individual spaces for better block assembly (see BlockSparseMatrixAssembler). Takes three parameters:
NB: Number of assembly blocks
SB: Size of each assembly block, as a Tuple.
P : Permutation of the variables of the multifield space when assembling, as a Tuple.
The resulting MultiFieldFEFunction is in the space (in particular it fulfills Dirichlet BCs even in the case that the given cell field does not fulfill them)
like interpolate, but also compute new degrees of freedom for the dirichlet component. The resulting MultiFieldFEFunction does not necessary belongs to the underlying space
The resulting MultiFieldFEFunction is in the space (in particular it fulfills Dirichlet BCs even in the case that the given cell field does not fulfill them)
like interpolate, but also compute new degrees of freedom for the dirichlet component. The resulting MultiFieldFEFunction does not necessary belongs to the underlying space
where $\overline{\alpha_{M}} = 1 - \alpha_{M}$, $\overline{\alpha_{F}} = 1 - \alpha_{F}$, $\overline{\gamma} = 1 - \gamma$ and $\overline{\beta} = \frac{1}{2}(1 - 2 \beta)$. Here again, we immediately see that $\boldsymbol{u}_{n+1}$ satisfies the recurrence
These conditions are hard to examine analytically, but one can verify that this scheme is at least of order $1$. Second-order is achieved by setting $\gamma = \frac{1}{2} - \alpha_{M} + \alpha_{F}$.
It is easier to consider the limit cases $|z| \to 0$ and $|z| \to +\infty$ and look at the eigenvalues of the amplification matrix.
When $|z| \to 0$, we find $\rho(\boldsymbol{A}(z)) = \max\left\{1, \left|\frac{\alpha_{M}}{1 - \alpha_{M}}\right|\right\}$.
For all these eigenvalues to have a modulus smaller than one, we need $\alpha_{M} \leq \frac{1}{2}$, $\alpha_{F} \leq \frac{1}{2}$, $\gamma \geq \frac{1}{2}$, i.e. $\alpha_{F} \geq \alpha_{M}$ and $\beta \geq \frac{1}{2} \gamma$. Since dissipation of high-frequency is maximised when the eigenvalues are real at infinity, we also impose $\beta = \frac{1}{16} (1 + 2 \gamma)^{2}$, i.e. $\beta = \frac{1}{4} (1 - \alpha_{M} + \alpha_{F})^{2}$.
This method was also designed to damp high-frequency perturbations so it is common practice to parameter this scheme in terms of its spectral radius.
The Hilbert-Huges-Taylor-$\alpha$ (HHT-$\alpha$) method is obtained by setting $\alpha_{M} = 0$, $\alpha_{F} = \frac{1 - \rho_{\infty}}{1 + \rho_{\infty}}$.
The Wood-Bossak-Zienkiewicz-$\alpha$ (WBZ-$\alpha$) method is recovered by setting $\alpha_{F} = 0$ and $\alpha_{M} = \frac{\rho_{\infty} - 1}{\rho_{\infty} + 1}$.
The standard generalised-$\alpha$ method is obtained by setting $\alpha_{M} = \frac{2 \rho_{\infty - 1}}{\rho_{\infty} + 1}$, $\alpha_{F} = \frac{\rho_{\infty}}{\rho_{\infty} + 1}$.
The Newmark method corresponds to $\alpha_{F} = \alpha_{M} = 0$. In this case, the values of $\beta$ and $\gamma$ are usually chosen as $\beta = 0$, $\gamma = \frac{1}{2}$ (explicit central difference scheme), or $\beta = \frac{1}{4}$ and $\gamma = \frac{1}{2}$ (midpoint rule).
where $\overline{\alpha_{M}} = 1 - \alpha_{M}$, $\overline{\alpha_{F}} = 1 - \alpha_{F}$, $\overline{\gamma} = 1 - \gamma$ and $\overline{\beta} = \frac{1}{2}(1 - 2 \beta)$. Here again, we immediately see that $\boldsymbol{u}_{n+1}$ satisfies the recurrence
These conditions are hard to examine analytically, but one can verify that this scheme is at least of order $1$. Second-order is achieved by setting $\gamma = \frac{1}{2} - \alpha_{M} + \alpha_{F}$.
It is easier to consider the limit cases $|z| \to 0$ and $|z| \to +\infty$ and look at the eigenvalues of the amplification matrix.
When $|z| \to 0$, we find $\rho(\boldsymbol{A}(z)) = \max\left\{1, \left|\frac{\alpha_{M}}{1 - \alpha_{M}}\right|\right\}$.
For all these eigenvalues to have a modulus smaller than one, we need $\alpha_{M} \leq \frac{1}{2}$, $\alpha_{F} \leq \frac{1}{2}$, $\gamma \geq \frac{1}{2}$, i.e. $\alpha_{F} \geq \alpha_{M}$ and $\beta \geq \frac{1}{2} \gamma$. Since dissipation of high-frequency is maximised when the eigenvalues are real at infinity, we also impose $\beta = \frac{1}{16} (1 + 2 \gamma)^{2}$, i.e. $\beta = \frac{1}{4} (1 - \alpha_{M} + \alpha_{F})^{2}$.
This method was also designed to damp high-frequency perturbations so it is common practice to parameter this scheme in terms of its spectral radius.
The Hilbert-Huges-Taylor-$\alpha$ (HHT-$\alpha$) method is obtained by setting $\alpha_{M} = 0$, $\alpha_{F} = \frac{1 - \rho_{\infty}}{1 + \rho_{\infty}}$.
The Wood-Bossak-Zienkiewicz-$\alpha$ (WBZ-$\alpha$) method is recovered by setting $\alpha_{F} = 0$ and $\alpha_{M} = \frac{\rho_{\infty} - 1}{\rho_{\infty} + 1}$.
The standard generalised-$\alpha$ method is obtained by setting $\alpha_{M} = \frac{2 \rho_{\infty - 1}}{\rho_{\infty} + 1}$, $\alpha_{F} = \frac{\rho_{\infty}}{\rho_{\infty} + 1}$.
The Newmark method corresponds to $\alpha_{F} = \alpha_{M} = 0$. In this case, the values of $\beta$ and $\gamma$ are usually chosen as $\beta = 0$, $\gamma = \frac{1}{2}$ (explicit central difference scheme), or $\beta = \frac{1}{4}$ and $\gamma = \frac{1}{2}$ (midpoint rule).
The implicit operator defined by the implicit residual is considered stiff and is meant to be solved implicitly,
The explicit operator defined by the explicit residual is considered non-stiff and is meant to be solved explicitly.
Important
The explicit operator must have one order less than the implicit operator, so that the mass term of the global operator is fully contained in the implicit operator.
The implicit operator defined by the implicit residual is considered stiff and is meant to be solved implicitly,
The explicit operator defined by the explicit residual is considered non-stiff and is meant to be solved explicitly.
Important
The explicit operator must have one order less than the implicit operator, so that the mass term of the global operator is fully contained in the implicit operator.
Linear stage operator representing res(x) = J(t, us) x + r(t, us) = 0, where x is the stage unknown and us denotes the point where the residual of the ODE is to be evaluated.
struct NonlinearStageOperator <: StageOperator end
Nonlinear stage operator representing res(x) = residual(t, us(x)...) = 0, where x is the stage unknown and us(x) denotes the point where the residual of the ODE is to be evaluated. It is assumed that the coordinates of us(x) are linear in x, and the coefficients in front of x called ws are scalar, i.e. ws[k] = d/dx us[k](x) is a scalar constant.
Wrapper that transforms a TransientFEOperator into an ODEOperator, i.e. takes residual(t, uh, ∂t[uh], ..., ∂t^N[uh], vh) and returns residual(t, us), where us[k] = ∂t^k[us] and uf represents the free values of uh.
Wrapper around an ODEOperator and ODESolver that represents the solution at a set of time steps. It is an iterator that computes the solution at each time step in a lazy fashion when accessing the solution.
An ODESolver is a map that update state vectors. These state vectors are created at the first iteration from the initial conditions, and are then converted back into the evaluation of the solution at the current time step.
In the simplest case, the state vectors correspond to the first N-1 time derivatives of u at time t_n, where N is the order of the ODEOperator, but some solvers rely on other state variables (values at previous times, higher-order derivatives...).
Linear stage operator representing res(x) = J(t, us) x + r(t, us) = 0, where x is the stage unknown and us denotes the point where the residual of the ODE is to be evaluated.
struct NonlinearStageOperator <: StageOperator end
Nonlinear stage operator representing res(x) = residual(t, us(x)...) = 0, where x is the stage unknown and us(x) denotes the point where the residual of the ODE is to be evaluated. It is assumed that the coordinates of us(x) are linear in x, and the coefficients in front of x called ws are scalar, i.e. ws[k] = d/dx us[k](x) is a scalar constant.
Wrapper that transforms a TransientFEOperator into an ODEOperator, i.e. takes residual(t, uh, ∂t[uh], ..., ∂t^N[uh], vh) and returns residual(t, us), where us[k] = ∂t^k[us] and uf represents the free values of uh.
Wrapper around an ODEOperator and ODESolver that represents the solution at a set of time steps. It is an iterator that computes the solution at each time step in a lazy fashion when accessing the solution.
An ODESolver is a map that update state vectors. These state vectors are created at the first iteration from the initial conditions, and are then converted back into the evaluation of the solution at the current time step.
In the simplest case, the state vectors correspond to the first N-1 time derivatives of u at time t_n, where N is the order of the ODEOperator, but some solvers rely on other state variables (values at previous times, higher-order derivatives...).
TimeSpaceFunction allows for convenient ways to apply differential operators to functions that depend on time and space. More precisely, if f is a function that, to a given time, returns a function of space (i.e. f is evaluated at time t and position x via f(t)(x)), then F = TimeSpaceFunction(f) supports the following syntax:
op(F): a TimeSpaceFunction representing both t -> x -> op(f)(t)(x) and (t, x) -> op(f)(t)(x),
op(F)(t): a function of space representing x -> op(f)(t)(x)
op(F)(t, x): the quantity op(f)(t)(x) (this notation is equivalent to op(F)(t)(x)),
for all spatial and temporal differential operator, i.e. op in (time_derivative, gradient, symmetric_gradient, divergence, curl, laplacian) and their symbolic aliases (∂t, ∂tt, ∇, ...).
abstract type TransientFEOperator <: GridapType end
Transient version of FEOperator corresponding to a residual of the form
residual(t, u, v) = 0,
where residual is linear in v. Time derivatives of u can be included by using the ∂t operator.
Important
For now, the residual and jacobians cannot be directly computed on a TransientFEOperator. They have to be evaluated on the corresponding algebraic operator, which is an ODEOperator. As such, TransientFEOperator is not exactly a subtype of FEOperator, but rather at the intersection of FEOperator and ODEOperator. This is because the ODEOperator works with vectors and it is optimised to take advantage of constant forms.
abstract type TransientFESolution <: GridapType end
Wrapper around a TransientFEOperator and ODESolver that represents the solution at a set of time steps. It is an iterator that computes the solution at each time step in a lazy fashion when accessing the solution.
abstract type TransientIMEXFEOperator <: TransientFEOperator end
Implicit-Explicit decomposition of a residual defining a TransientFEOperator:
residual(t, u, v) = implicit_residual(t, u, v)
- + explicit_residual(t, u, v),
where
The implicit operator defined by the implicit residual is considered stiff and is meant to be solved implicitly,
The explicit operator defined by the explicit residual is considered non-stiff and is meant to be solved explicitly.
Both the implicit and explicit residuals are linear in v.
Important
The explicit operator must have one order less than the implicit operator, so that the mass term of the global operator is fully contained in the implicit operator.
struct TransientLinearFEOpFromWeakForm <: TransientFEOperator end
Transient FEOperator defined by a transient weak form
residual(t, u, v) = ∑_{0 ≤ k ≤ N} form_k(t, ∂t^k[u], v) - res(t, v) = 0,
where N is the order of the operator, form_k is linear in ∂t^k[u] and does not depend on the other time derivatives of u, and the form_k and res are linear in v.
For convenience, the form corresponding to order k has to be written as a function of ∂t^k[u], i.e. as a linear form, and the residual as a function of t and v only.
TimeSpaceFunction allows for convenient ways to apply differential operators to functions that depend on time and space. More precisely, if f is a function that, to a given time, returns a function of space (i.e. f is evaluated at time t and position x via f(t)(x)), then F = TimeSpaceFunction(f) supports the following syntax:
op(F): a TimeSpaceFunction representing both t -> x -> op(f)(t)(x) and (t, x) -> op(f)(t)(x),
op(F)(t): a function of space representing x -> op(f)(t)(x)
op(F)(t, x): the quantity op(f)(t)(x) (this notation is equivalent to op(F)(t)(x)),
for all spatial and temporal differential operator, i.e. op in (time_derivative, gradient, symmetric_gradient, divergence, curl, laplacian) and their symbolic aliases (∂t, ∂tt, ∇, ...).
abstract type TransientFEOperator <: GridapType end
Transient version of FEOperator corresponding to a residual of the form
residual(t, u, v) = 0,
where residual is linear in v. Time derivatives of u can be included by using the ∂t operator.
Important
For now, the residual and jacobians cannot be directly computed on a TransientFEOperator. They have to be evaluated on the corresponding algebraic operator, which is an ODEOperator. As such, TransientFEOperator is not exactly a subtype of FEOperator, but rather at the intersection of FEOperator and ODEOperator. This is because the ODEOperator works with vectors and it is optimised to take advantage of constant forms.
abstract type TransientFESolution <: GridapType end
Wrapper around a TransientFEOperator and ODESolver that represents the solution at a set of time steps. It is an iterator that computes the solution at each time step in a lazy fashion when accessing the solution.
abstract type TransientIMEXFEOperator <: TransientFEOperator end
Implicit-Explicit decomposition of a residual defining a TransientFEOperator:
residual(t, u, v) = implicit_residual(t, u, v)
+ + explicit_residual(t, u, v),
where
The implicit operator defined by the implicit residual is considered stiff and is meant to be solved implicitly,
The explicit operator defined by the explicit residual is considered non-stiff and is meant to be solved explicitly.
Both the implicit and explicit residuals are linear in v.
Important
The explicit operator must have one order less than the implicit operator, so that the mass term of the global operator is fully contained in the implicit operator.
struct TransientLinearFEOpFromWeakForm <: TransientFEOperator end
Transient FEOperator defined by a transient weak form
residual(t, u, v) = ∑_{0 ≤ k ≤ N} form_k(t, ∂t^k[u], v) - res(t, v) = 0,
where N is the order of the operator, form_k is linear in ∂t^k[u] and does not depend on the other time derivatives of u, and the form_k and res are linear in v.
For convenience, the form corresponding to order k has to be written as a function of ∂t^k[u], i.e. as a linear form, and the residual as a function of t and v only.
Compute the jacobian of the residual of the ODEOperator with respect to all time derivatives, weighted by some factors ws.
The weights are ordered by increasing order of time derivative, i.e. the first weight corresponds to ∂residual / ∂u and the last to ∂residual / ∂(d^N u / dt^N).
Compute the jacobian of the residual of the ODEOperator with respect to all time derivatives, weighted by some factors ws.
The weights are ordered by increasing order of time derivative, i.e. the first weight corresponds to ∂residual / ∂u and the last to ∂residual / ∂(d^N u / dt^N).
Allocate a jacobian matrix for the ODEOperator and compute the jacobian of the residual of the ODEOperator with respect to all time derivatives, weighted by some factors ws.
The weights are ordered by increasing order of time derivative, i.e. the first weight corresponds to ∂residual / ∂u and the last to ∂residual / ∂(d^N u / dt^N).
Allocate a jacobian matrix for the ODEOperator and compute the jacobian of the residual of the ODEOperator with respect to all time derivatives, weighted by some factors ws.
The weights are ordered by increasing order of time derivative, i.e. the first weight corresponds to ∂residual / ∂u and the last to ∂residual / ∂(d^N u / dt^N).
Return the ODEOperatorType of the operator defined by an IMEX decomposition. This function should be called in the constructors of concrete IMEX operators.
Return the ODEOperatorType of the operator defined by an IMEX decomposition. This function should be called in the constructors of concrete IMEX operators.
Check whether two operators can make a valid IMEX operator decomposition. This function should be called in the constructors of concrete IMEX operators.
Check whether two operators can make a valid IMEX operator decomposition. This function should be called in the constructors of concrete IMEX operators.
Add the jacobian of the residual of the ODEOperator with respect to all time derivatives, weighted by some factors ws.
The weights are ordered by increasing order of time derivative, i.e. the first weight corresponds to ∂residual / ∂u and the last to ∂residual / ∂(d^N u / dt^N).
Add the jacobian of the residual of the ODEOperator with respect to all time derivatives, weighted by some factors ws.
The weights are ordered by increasing order of time derivative, i.e. the first weight corresponds to ∂residual / ∂u and the last to ∂residual / ∂(d^N u / dt^N).