refinement

This submodule’s purpose is to provide support functions for mesh refinement, namely the ones involved in post-solver refinement. Currently heavily work in progress.

Utilities

For computing Legendre decomposition of functions for error estimation, the function compute_legendre_coefficients() is provided. It is called by compute_legendre_error_estimates(), which is used to estimate the \(L^2\) error estimate and the approximate increase in the \(L^2\) error due to order drop with \(h\)-refinement.

mfv2d.refinement.compute_legendre_error_estimates(order_1: int, order_2: int, xi: npt.NDArray[np.float64], eta: npt.NDArray[np.float64], w: npt.NDArray[np.float64], det: npt.NDArray[np.float64], u: npt.NDArray[np.float64], err: npt.NDArray[np.float64]) → tuple[float, float]

Compute error estimates based on Legendre coefficients.

Parameters:

• order_1 (int) – Order of the element in the first direction.

• order_2 (int) – Order of the element in the second direction.

• xi (array) – Positions of integration nodes in the first direction.

• eta (array) – Positions of integration nodes in the second direction.

• w (array) – Values of integration weights.

• w – Values of the Jacobian determinant.

• u (array) – Solution values at the integration nodes.

• err (array) – Error estimate values at the integration nodes.

Returns:

• l2_error_square (float) – Square of the error estimate in the L2 norm.

• h_ref_cost (float) – Estimate of the cost for h-refinemnt.

mfv2d.refinement.compute_legendre_coefficients(order_1: int, order_2: int, nodes_xi: npt.NDArray[np.float64], nodes_eta: npt.NDArray[np.float64], weighted_function: npt.NDArray[np.float64], det: npt.NDArray[np.float64]) → npt.NDArray[np.float64]

Compute Legendre coefficients from function values at integration nodes.

Parameters:

• order_1 (int) – Order of the coefficients in the first direction.

• order_2 (int) – Order of the coefficients in the second direction.

• weighted_function (array) – Product of function and integration weight for all integration points multiplied by the Jacobian determinant.

• det (array) – Determinant of the Jacobian.

Returns:

Array of coefficients for Legendre basis.

Return type:

(order_1 + 1, order_2 + 1) array

Refinement Settings

Refinement settings are specified using the type RefinementSettings. The two most important members of it are the error_estimate and the refinement_limit members.

class mfv2d.refinement.RefinementSettings(error_estimate: ErrorEstimateCustom | ErrorEstimateLocalInverse | ErrorEstimateL2OrderReduction | ErrorEstimateExplicit | ErrorEstimateVMS, refinement_limit: RefinementLimitUnknownCount | RefinementLimitElementCount | RefinementLimitErrorValue, h_refinement_ratio: float = 0.0, report_error_distribution: bool = False, report_order_distribution: bool = False, upper_order_limit: int | None = None, lower_order_limit: int | None = None)

Settings pertaining to refinement of a mesh.

error_estimate: ErrorEstimateCustom | ErrorEstimateLocalInverse | ErrorEstimateL2OrderReduction | ErrorEstimateExplicit | ErrorEstimateVMS

How the error ought to be estimated.

h_refinement_ratio: float = 0.0

Ratio between element error and h-refinement cost where refinement can happen.

lower_order_limit: int | None = None

Elements bellow this order can not be h-refined.

refinement_limit: RefinementLimitUnknownCount | RefinementLimitElementCount | RefinementLimitErrorValue

Limit for mesh refinement.

report_error_distribution: bool = False

Should the error distribution be reported.

report_order_distribution: bool = False

Should the order distribution be reported.

upper_order_limit: int | None = None

Elements which reach this can only be h-refined.

Refinement Limits

To specify when the refineemnt should stop, several different dataclasses are available. RefinementLimitUnknownCount limits it based on number of new degrees of freedom that are introduced, RefinementLimitElementCount is based on the number of elements refined, added, and RefinementLimitErrorValue is based on the error value of elements.

class mfv2d.refinement.RefinementLimitUnknownCount(maximum_fraction: float, maximum_count: int)

Limit refinement based on change in number of degrees of freedom.

maximum_count: int

maximum_fraction: float

class mfv2d.refinement.RefinementLimitElementCount(maximum_fraction: float, maximum_count: int)

Limit refinement based on the number of elements.

maximum_count: int

maximum_fraction: float

class mfv2d.refinement.RefinementLimitErrorValue(minimum_fraction: float, minimum_value: float)

Limit refinement based on value of the errror.

minimum_fraction: float

minimum_value: float

Error Estimates

Error can be specified in different ways. These differ in time and accuracy.

class mfv2d.refinement.ErrorEstimateCustom(required_forms: Sequence[KFormUnknown], error_calculation_function: ErrorCalculationFunctionFull | ErrorCalculationFunctionSimple, reconstruction_orders: tuple[int, int] | None = None)

Settings for refinement that is to be done with user-defined function.

error_calculation_function: ErrorCalculationFunctionFull | ErrorCalculationFunctionSimple

Function called to calculate error estimate and h-refinement cost.

reconstruction_orders: tuple[int, int] | None = None

Order at which error should be reconstructed.

required_forms: Sequence[KFormUnknown]

Forms that are needed by the error calculation function.

class mfv2d.refinement.ErrorEstimateLocalInverse(target_form: KFormUnknown, order_increase: int, strong_forms: Sequence[KFormUnknown] = ())

Settings for refinement that is based on local inverse.

order_increase: int

Order at which residual should be reconstructed prior to inversion.

strong_forms: Sequence[KFormUnknown] = ()

Forms for which the boundary conditions on each element must be given strongly.

target_form: KFormUnknown

Error of this form is used as a guid for refinement.

class mfv2d.refinement.ErrorEstimateL2OrderReduction(target_form: KFormUnknown, order_drop: int, alternative: Literal['ignore', 'prioritize'] = 'prioritize')

Reduces order of solution in order to estimate error.

alternative: Literal['ignore', 'prioritize'] = 'prioritize'

Alternative measure to use whenever the current measure can not be applied due to insufficient element order.

• “ignore” means that elements with order less than order_drop will be skipped,

• “prioritize” meanst that elements with order less than order_drop will receive an infinite error estimate.

order_drop: int

Order of an element is dropped by this much.

target_form: KFormUnknown

Error of this form is used as a guid for refinement.

class mfv2d.refinement.ErrorEstimateExplicit(target_form: KFormUnknown, solution_estimate: Function2D, reconstruction_orders: tuple[int, int] | None = None)

Estimate error based on an explicit function.

reconstruction_orders: tuple[int, int] | None = None

Order at which error should be reconstructed.

solution_estimate: Function2D

Function of poistion, which is used to estimate the exact solution from.

target_form: KFormUnknown

Form for which the error is being computed.

class mfv2d.refinement.ErrorEstimateVMS(target_form: KFormUnknown, symmetric_system: KFormSystem, nonsymmetric_system: KFormSystem, order_increase: int, max_iters: int, atol: float, rtol: float)

Estimate error based on variational multi-scale approach.

atol: float

When the largest update in the fine scales drops bellow this value, consider it converged.

max_iters: int

Number of iterations to improve the guess of fine scales.

nonsymmetric_system: KFormSystem

Non-symmetric system based on which the fine-scale function is computed.

order_increase: int

Order increase for the Green’s function and residual.

rtol: float

When the largest update in the fine scalse drops bellow the largest value of fine scale degrees of freedom scaled by this value, consider it converged.

symmetric_system: KFormSystem

Symmetric system based on which the Green’s function is computed.

target_form: KFormUnknown

Form for which the error is being computed.

Each of these has its respective error estimation function:

• ErrorEstimateCustom has error_estimate_with_custom_estimator()

• ErrorEstimateLocalInverse has error_estimate_with_local_inversion()

• ErrorEstimateL2OrderReduction has error_estimate_with_order_reduction()

• ErrorEstimateExplicit has error_estimate_with_explicit_solution()

• ErrorEstimateVMS has error_estimate_with_vms()

mfv2d.refinement.error_estimate_with_custom_estimator(leaf_count: int, solution: npt.NDArray[np.float64], element_offsets: npt.NDArray[np.uint32], required_unknowns: Sequence[KFormUnknown], form_specs: ElementFormSpecification, error_calculation_function: CustomErrorFunction, element_fem_spaces: Sequence[ElementFemSpace2D], recon_order_1: int | None, recon_order_2: int | None) → tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]

Compute element error estimates using a custom, user-supplied function.

Parameters:

• leaf_count (int) – Number of leaves in the mesh.

• solution (array) – Full solution vector.

• element_offsets (array) – Array of offsets to the beginning of each element for the solution vector.

• dof_offsets (array) – Offsets of each unknown within each element.

• required_unknowns (Sequence of KFormUnknown) – Unknowns which are needed by the user function.

• unknown_forms (Sequence of KFormUnknown) – Unknown forms in the order they appear in the system.

• error_calculation_function (CustomErrorFunction) – Error function that is used to compute the error estimates.

• elemenebt_fem_spaces (Sequence of ElementFemSpace2D) – Element FEM spaces.

• recon_order_1 (int or None) – Order at which all element should be reconstructed. If not present, the element integration order is used.

recon_order_2int or None

Order at which all element should be reconstructed. If not present, the element integration order is used.

Returns:

• element_error (array) – Array with error estimates for each element.

• href_cost (array) – Array with estimates of increase in error due to h-refinement.

mfv2d.refinement.error_estimate_with_local_inversion(mesh: Mesh, solution: npt.NDArray[np.float64], element_offsets: npt.NDArray[np.uint32], boundary_conditions: Sequence[BoundaryCondition2DSteady], element_fem_spaces: Sequence[ElementFemSpace2D], recon_order: int, basis_cache: FemCache, system: KFormSystem, compiled: CompiledSystem, unknown_target: KFormUnknown, strongly_zeroed: Sequence[KFormUnknown], constrained: Sequence[KFormUnknown]) → tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]

Compute element error estimates using a local inversion.

Returns:

• element_error (array) – Array with error estimates for each element.

• href_cost (array) – Array with estimates of increase in error due to h-refinement.

mfv2d.refinement.error_estimate_with_order_reduction(solution: npt.NDArray[np.float64], element_offsets: npt.NDArray[np.uint32], element_fem_spaces: Sequence[ElementFemSpace2D], reduction_order: int, basis_cache: FemCache, unknown_forms: ElementFormSpecification, unknown_target: KFormUnknown, alternative: Literal['ignore', 'prioritize']) → tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]

Compute element error estimates using a local inversion.

Returns:

• element_error (array) – Array with error estimates for each element.

• href_cost (array) – Array with estimates of increase in error due to h-refinement.

mfv2d.refinement.error_estimate_with_explicit_solution(solution: npt.NDArray[np.float64], element_offsets: npt.NDArray[np.uint32], required_unknown: KFormUnknown, form_specs: ElementFormSpecification, soulution_calculation_function: Function2D, element_fem_spaces: Sequence[ElementFemSpace2D], recon_order_1: int | None, recon_order_2: int | None, basis_cache: FemCache) → tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]

Compute element error estimates using a custom, user-supplied function.

Parameters:

• leaf_count (int) – Number of leaves in the mesh.

• solution (array) – Full solution vector.

• element_offsets (array) – Array of offsets to the beginning of each element for the solution vector.

• dof_offsets (array) – Offsets of each unknown within each element.

• required_unknowns (Sequence of KFormUnknown) – Unknowns which are needed by the user function.

• unknown_forms (Sequence of KFormUnknown) – Unknown forms in the order they appear in the system.

• error_calculation_function (CustomErrorFunction) – Error function that is used to compute the error estimates.

• elemenebt_fem_spaces (Sequence of ElementFemSpace2D) – Element FEM spaces.

• recon_order_1 (int or None) – Order at which all element should be reconstructed. If not present, the element integration order is used.

recon_order_2int or None

Order at which all element should be reconstructed. If not present, the element integration order is used.

Returns:

• element_error (array) – Array with error estimates for each element.

• href_cost (array) – Array with estimates of increase in error due to h-refinement.

mfv2d.refinement.error_estimate_with_vms(mesh: Mesh, leaf_indices: Sequence[int], solution: npt.NDArray[np.float64], element_offsets: npt.NDArray[np.uint32], boundary_conditions: Sequence[BoundaryCondition2DSteady], element_fem_spaces: Sequence[ElementFemSpace2D], recon_order: int, basis_cache: FemCache, system: KFormSystem, compiled: CompiledSystem, symmetric: KFormSystem, nonsymmetric: KFormSystem, unknown_target: KFormUnknown, constrained_forms: Sequence[tuple[float, KFormUnknown]], atol: float, rtol: float, max_iters: int) → tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]

Compute element error estimates using a local inversion.

Returns:

• element_error (array) – Array with error estimates for each element.

• href_cost (array) – Array with estimates of increase in error due to h-refinement.

Doing the Refinement

After the solver obtains a solution for the problem, refinement is performed by invoking the function perform_mesh_refinement(). This takes values computed by the solver, as well as all that is specified through RefinementSettings. After it calls the error estimation function of the error estimate specified, refine_mesh_based_on_error() is called to perform the refinement based on the computed error measures.

mfv2d.refinement.perform_mesh_refinement(mesh: Mesh, solution: npt.NDArray[np.float64], element_offsets: npt.NDArray[np.uint32], system: KFormSystem, error_estimator: ErrorEstimate, h_refinement_ratio: float, refinement_limit: RefinementLimit, report_error_distribution: bool, element_fem_spaces: Sequence[ElementFemSpace2D], boundary_conditions: Sequence[BoundaryCondition2DSteady], basis_cache: FemCache, order_limit: int | None, lower_order_limit: int | None, constrained: Sequence[tuple[float, KFormUnknown]]) → tuple[Mesh, npt.NDArray[np.float64], npt.NDArray[np.float64]]

Perform a round of mesh refinement.

Parameters:

• mesh (Mesh) – Mesh on which to perform refinement on.

• solution (array) – Solution for degrees of freedom.

• element_offsets (array) – Array of offsets of degrees of freedom.

• system (KFormSystem) – System that is being solved.

• error_estimator (ErrorEstimate) – Error estimator for the system.

• h_refinement_ratio (float) – Ratio of h-refinement error estimate to total element error estimate, at which h-refinement should occur.

• refinemenet_limit (RefinementLimit) – When should the refinement stop.

• report_error_distribution (bool) – Should the error distribution be reported to the terminal.

• element_fem_spaces (Sequence of ElementFemSpace2D) – FEM spaces of the elements.

• boundary_conditions (Sequence of BoundaryCondition2DSteady) – Sequence of boundary conditions.

• basis_cache (FemCache) – Cache of basis.

• leaf_order_mapping (dict of int -> int) – Dictionary that maps leaf element indices to their order in the system.

Returns:

• Mesh – Newly created refined mesh.

• array – Array of error estimates for each element.

• array – Array of h-refinement error estimates for each element.

mfv2d.refinement.refine_mesh_based_on_error(mesh: Mesh, total_unknowns: int, h_refinement_ratio: float, refinement_limit: RefinementLimit, form_specs: ElementFormSpecification, leaf_indices: npt.NDArray[np.uint32], element_error: npt.NDArray[np.float64], href_cost: npt.NDArray[np.float64], order_limit: int | None, lower_order_limit: int | None) → Mesh

Refine the given mesh based on given element error and h-refinement cost.