Python API

The problem

class dualmesh.Problem(mesh, method='dmcdm', coordinates='cartesian', boundary_gradient='first_order')[source]

A boundary value problem discretized on a primal mesh and its dual mesh.

Parameters:

  • mesh – The primal mesh of finite elements.

  • method (str) – The discretization. "dmcdm" is the dual mesh control domain method (the default); "fem" is the Galerkin finite element method; "hfvm" is the vertex-centred finite volume method, which uses the same control domains as the dual mesh method but two-point gradients at their interfaces (the half-control volume formulation of Reddy, Chapter 3); "zfvm" is the cell-centred finite volume method, with one unknown per element and one per boundary face (the zero-thickness control volume formulation, and the layout used by OpenFOAM). The rest of the problem definition is identical for all four, which makes them directly comparable.

  • boundary_gradient (str) – Only for "zfvm": "first_order" (the two-point difference between the cell and the boundary node, the default) or "second_order" (the one-sided quadratic of Eq. (3.2.14) of the book, through the boundary node and the two nearest cells).

  • coordinates (str) – "cartesian", "axisymmetric" (the integrals carry the factor \(2 \pi r\), with \(r\) the first coordinate), or "spherical" (factor \(4 \pi r^2\), one-dimensional meshes).

set_num_threads(num_threads)[source]

Number of threads used by the assembly loops.

0 means the OpenMP default, which is the value of the OMP_NUM_THREADS environment variable or, if that is unset, the number of cores. The setting is a request: the assembly falls back to one thread when the library was built without OpenMP, when the mesh is too small for threading to pay for itself, or when any kernel, boundary condition, material or function of the problem is defined in Python, because calling back into the interpreter needs the global interpreter lock. effective_threads() reports the number actually used.

Parameters:

num_threads (int)

Return type:

None

effective_threads()[source]

The number of threads the assembly will actually use.

Return type:

int

property thread_safe: bool

False when some object of the problem is defined in Python, which forces the assembly onto one thread.

property is_cell_centered: bool

True when the unknowns sit at cell centroids instead of mesh nodes.

entity_points()[source]

Positions of the degrees of freedom, one row each, in the order of values().

For the dual mesh, finite element and vertex-centred finite volume methods these are the mesh nodes. For the cell-centred finite volume method they are the cell centroids first and then the boundary face centroids.

Return type:

numpy.ndarray

boundary_entities(boundary)[source]

Degree of freedom indices that carry the values of a boundary.

Parameters:

boundary (str)

Return type:

list[int]

add_variable(name, blocks=(), initial_condition=None)[source]

Add a nodal unknown; returns its index.

Parameters:
Return type:

int

add_function(name, function)[source]

Register a named function of (x, y, z, t) (or a constant).

Parameters:

name (str)

Return type:

None

add_kernel(kernel, name=None, **parameters)[source]

Add a kernel by registered type name, or an object built in Python.

Parameters:

name (str | None)

add_boundary_condition(condition, name=None, **parameters)[source]

Add a boundary condition (essential or natural).

Parameters:

name (str | None)

add_material(material, name=None, **parameters)[source]

Add a material (a provider of named properties).

Parameters:

name (str | None)

add_point_source(source='PointSource', name=None, **parameters)[source]

Add a concentrated nodal source (point force or point heat source).

Parameters:

name (str | None)

add_nodal_load(source='PointSource', name=None, **parameters)

Add a concentrated nodal source (point force or point heat source).

Parameters:

name (str | None)

initialize()[source]

Resolve all objects (called automatically by solve()).

Return type:

None

solve(**options)[source]

Solve the steady problem.

Keyword arguments are solver options: nonlinear_solver ("newton", "picard", or "linear"), max_iterations, relative_tolerance, absolute_tolerance, step_tolerance, relaxation (the acceleration parameter of direct iteration), load_factors (load stepping), linear_solver ("automatic", the default, "lu", "bicgstab", "gmres" or "cg"), preconditioner ("ilu", the default, "ilut", "jacobi" or "none"), linear_tolerance, linear_max_iterations, gmres_restart, verbose, and error_on_divergence.

"automatic" factorises the system directly where that is cheap, which is always in one dimension, up to \(10^5\) unknowns in two and up to a few thousand in three, and otherwise uses BiCGSTAB preconditioned by an incomplete LU factorisation, falling back to the direct solver if the iteration does not converge. On a three-dimensional mesh the iteration is typically ten to fifty times faster than the direct solver, because the direct factors of a 3D problem fill in far more.

Return type:

dualmesh._core.SolveResult

solve_transient(end_time, dt, start_time=0.0, theta=1.0, output_interval=0, output_file_base='', time_stepper='fixed', dt_min=0.0, dt_max=0.0, growth_factor=2.0, cutback_factor=0.5, error_tolerance=0.001, optimal_iterations=4, iteration_window=2, max_rejected_steps=10, **options)[source]

March the solution forward in time with the theta method.

The theta method weights the steady part of the residual between the old state and the new one,

\[R_{\text{time}}(U^{n+1}) + \theta R_{\text{steady}}(U^{n+1}) + (1 - \theta) R_{\text{steady}}(U^{n}) = 0,\]

so theta=1 is backward Euler, which is unconditionally stable and first-order accurate; theta=0.5 is the Crank-Nicolson, or midpoint, rule, which is unconditionally stable and second-order accurate but can ring on a sharp transient; and theta=0 is forward Euler, which is explicit in the steady terms and is stable only below a critical step.

Parameters:

  • end_time (float) – The interval to cover and the step to take. With an adaptive stepper dt is the first step rather than every step.

  • dt (float) – The interval to cover and the step to take. With an adaptive stepper dt is the first step rather than every step.

  • start_time (float) – The interval to cover and the step to take. With an adaptive stepper dt is the first step rather than every step.

  • theta (float) – The weight above.

  • output_interval (int) – Write a .vtu file every so many accepted steps.

  • output_file_base (str) – Write a .vtu file every so many accepted steps.

  • time_stepper (str) –

    "fixed" keeps the step, shortening only the last one so that the run lands exactly on end_time.

    "error" chooses the step from an estimate of the local truncation error. Each interval is advanced twice, once with one step and once with two half steps, and the difference between the two answers estimates the error of the coarse one by Richardson extrapolation. A step whose relative error exceeds error_tolerance is discarded and retried with a smaller step; an accepted step is followed by the largest step the estimate allows. The answer that is kept is the accurate one, from the two half steps. The estimate costs three nonlinear solves per accepted step, so use this when accuracy in time is what matters.

    "iteration" chooses the step from how hard the nonlinear solver worked. A step that converged in fewer than optimal_iterations - iteration_window iterations is followed by a larger one, a step that needed more than optimal_iterations + iteration_window by a smaller one, and a step that failed to converge is discarded and retried. It costs nothing beyond the solve and is the right choice when the difficulty is the nonlinearity rather than the accuracy.

  • dt_min (float) – Bounds on the step. A run that has to go below dt_min is reported as a failure rather than grinding to a halt. Zero means dt divided by one million, and the whole interval, respectively.

  • dt_max (float) – Bounds on the step. A run that has to go below dt_min is reported as a failure rather than grinding to a halt. Zero means dt divided by one million, and the whole interval, respectively.

  • growth_factor (float) – The most the step may grow between accepted steps, and the factor applied after a rejected one.

  • cutback_factor (float) – The most the step may grow between accepted steps, and the factor applied after a rejected one.

  • error_tolerance (float) – The target for the relative local error of one step.

  • optimal_iterations (int) – The iteration count the "iteration" stepper aims for, and the half-width of the band around it inside which the step is left alone.

  • iteration_window (int) – The iteration count the "iteration" stepper aims for, and the half-width of the band around it inside which the step is left alone.

  • max_rejected_steps (int) – How many times in a row a step may be rejected before the run is declared a failure.

  • **options – Passed to the nonlinear solver of every step; see solve().

Returns:

Besides the usual fields, time_steps counts the accepted steps, rejected_steps the discarded ones, and step_history holds the time reached and the step taken for each accepted step.

Return type:

SolveResult

set_time_step_callback(callback)[source]

Call callback(time, problem) after every converged time step.

Return type:

None

property mesh
property time: float
variable_index(name)[source]
Parameters:

name (str)

Return type:

int

values(variable)[source]

Values of a variable, one per degree of freedom entity.

For every method but "zfvm" these are nodal values indexed by node; for "zfvm" they are cell values followed by boundary face values. Use entity_points() for the matching coordinates.

Parameters:

variable (str)

Return type:

numpy.ndarray

set_values(variable, values)[source]
Parameters:

variable (str)

Return type:

None

solution()[source]

The full solution vector (node-major, variable-minor).

Return type:

numpy.ndarray

apply_initial_conditions()[source]
Return type:

None

reactions(variable, boundary)[source]

Secondary variables at the nodes of a boundary: [(node, value), ...].

For each node the value is the integral of the normal flux over the part of the boundary that belongs to that node’s control domain, which is the quantity Reddy denotes \(Q_I\) (a reaction, a heat flow, a force).

A node that lies on two boundaries (a corner) carries one reaction that covers its whole boundary portion, so it appears in both lists. When summing over several boundaries, collect the nodes first (dict(problem.reactions(...))) instead of adding the totals, or corner nodes are counted twice.

Parameters:

  • variable (str)

  • boundary (str)

total_reaction(variable, boundary)[source]

Sum of the secondary variables over a boundary.

Parameters:

  • variable (str)

  • boundary (str)

Return type:

float

sample(variable, points)[source]

Interpolate a variable at arbitrary points (NaN outside the mesh).

Parameters:

variable (str)

Return type:

numpy.ndarray

values_at_nodes(variable, nodes)[source]
Parameters:
Return type:

numpy.ndarray

nodes_where(predicate)[source]

Node indices whose coordinates satisfy predicate(x, y, z).

Return type:

list[int]

node_at(point, tolerance=1e-09)[source]

Index of the node closest to point (an error if none is within tolerance).

Parameters:

tolerance (float)

Return type:

int

gradient_at_centroids(variable)[source]
Parameters:

variable (str)

Return type:

numpy.ndarray

property_at_centroids(property_name)[source]
Parameters:

property_name (str)

Return type:

numpy.ndarray

kernel_flux_at_centroids(kernel_name)[source]
Parameters:

kernel_name (str)

Return type:

numpy.ndarray

error_indicator(variable)[source]

One error indicator per element, from gradient recovery.

The gradient of the computed solution jumps between elements. A smoother gradient is recovered by averaging the element gradients onto the nodes, weighted by the share of each element that belongs to the node’s control domain, and interpolating that nodal field back over the element. The indicator of an element is the square root of the integral over it of the squared difference between the two gradients. The recovered gradient is the more accurate of the two, so their difference measures the error in the computed one; this is the estimator of Zienkiewicz and Zhu (1987).

It is an indicator, not a bound. It says which elements carry most of the error, which is what mark_by_fraction() and its relatives need, and it does not certify the size of the error.

Parameters:

variable (str)

Return type:

numpy.ndarray

error_norms(variable, exact, exact_gradient=None, quadrature_points=0)[source]

The error of the computed field against a known exact solution.

Returns (l2, h1_seminorm): the \(L^2\) norm of \(u_h - u\) and the \(H^1\) seminorm, the \(L^2\) norm of \(\nabla u_h - \nabla u\). These are the norms in which the convergence theory of every method in the library is stated, so they are what a convergence study should measure.

exact is anything a parameter accepts: a number, an expression string, a ParsedFunction or a Python callable of (x, y, z, t). exact_gradient is a sequence of up to three of the same, one per component; missing components are taken as zero, and without it the seminorm is returned as nan.

\(u_h\) is the field the method actually represents: the element interpolation of the nodal values for dmcdm, fem and hfvm, and for zfvm the linear reconstruction \(U_c + G_c \\cdot (x - x_c)\) in every cell, whose gradient is the reconstructed cell gradient. The integrals use a Gauss rule of quadrature_points per direction on every element (zero chooses the polynomial order plus two, which over-integrates the leading term of the error) and include the coordinate factor. The sum runs on several threads unless one of the functions is a Python callable.

Parameters:

  • variable (str)

  • quadrature_points (int)

Return type:

tuple[float, float]

linear_system()[source]

The residual and the Jacobian of the steady problem at the current solution, as (residual, jacobian).

This is the system one Newton step of solve() solves, \(J\,\delta U = -R\): the prescribed boundary values are written into the solution first, and the rows (and columns) of the prescribed degrees of freedom are replaced by those of the identity. residual is a NumPy array and jacobian a SciPy compressed sparse column matrix, so the system can be handed to any solver or preconditioner that works with SciPy, for instance to compare linear solvers or to study the spectrum of a discretisation. It needs SciPy. The degree of freedom of variable v on entity i is i * num_variables + v.

integrate(variable)[source]
Parameters:

variable (str)

Return type:

float

boundary_flux_integral(kernel_name, boundary)[source]
Parameters:

  • kernel_name (str)

  • boundary (str)

Return type:

float

write_vtu(filename, cell_properties=())[source]

Write a VTK unstructured grid (readable by ParaView and VisIt).

Parameters:
Return type:

None

write_mesh_file(filename, file_format=None)[source]

Write the mesh and all nodal fields through meshio (Exodus, VTU, …).

Parameters:

  • filename (str)

  • file_format (str | None)

Return type:

None

write_csv(filename, variables=())[source]

Write nodal coordinates and values as comma-separated values.

Parameters:
Return type:

None

summary()[source]
Return type:

str

Meshes

class dualmesh.Mesh(*args, **kwargs)
dualmesh.generate_line_mesh(start=None, end=None, num_elements=None, bias=1.0, coordinates=None, element_type='Edge2')[source]

One-dimensional mesh of Edge2 or Edge3 elements.

Side sets and node sets "left" and "right" are created automatically. num_elements always counts elements, so an Edge3 mesh with num_elements=4 has four elements and nine nodes.

Parameters:
Return type:

dualmesh._core.Mesh

dualmesh.generate_rectangle_mesh(x_min=None, x_max=None, y_min=None, y_max=None, num_x_elements=None, num_y_elements=None, element_type='Quad4', diagonal='right', x_bias=1.0, y_bias=1.0, x_coordinates=None, y_coordinates=None)[source]

Structured two-dimensional mesh of Quad4, Tri3, Quad8, Quad9 or Tri6 elements. Quad8 works with method="fem" and method="zfvm" only.

Side sets "left", "right", "bottom" and "top" are created automatically. Pass x_coordinates/y_coordinates for a fully non-uniform grid, or the bounds together with the element counts and an optional geometric bias. The element counts always count elements, so asking for a quadratic type gives the same number of elements and more nodes, not fewer elements.

Parameters:
Return type:

dualmesh._core.Mesh

dualmesh.generate_box_mesh(x_min=None, x_max=None, y_min=None, y_max=None, z_min=None, z_max=None, num_x_elements=None, num_y_elements=None, num_z_elements=None, element_type='Hex8', x_bias=1.0, y_bias=1.0, z_bias=1.0, x_coordinates=None, y_coordinates=None, z_coordinates=None)[source]

Structured three-dimensional mesh of Hex8, Tet4, Wedge6, Pyramid5, Hex20, Hex27 or Tet10 elements. A Tet4 mesh has six tetrahedra per cell, a Wedge6 mesh two prisms split on the base diagonal, and a Pyramid5 mesh six pyramids meeting at an added centre node. Pyramid5 and Hex20 work with method="fem" and method="zfvm" only.

Side sets "left", "right", "bottom", "top", "back" and "front" are created automatically.

Parameters:
Return type:

dualmesh._core.Mesh

dualmesh.generate_annulus_mesh(inner_radius, outer_radius, num_radial_elements, num_angular_elements, start_angle=0.0, end_angle=90.0, element_type='Quad4', radial_bias=1.0, radial_coordinates=None)[source]

Mesh of an annular sector (angles in degrees, measured from the x axis).

Side sets: "inner" (r = inner_radius), "outer" (r = outer_radius), "start" (the start_angle edge) and "end" (the end_angle edge). This is the mesh used for thick pressurized cylinders and for the classical plate-with-a-hole problem.

Parameters:
Return type:

dualmesh._core.Mesh

dualmesh.read_mesh(filename, file_format=None, boundary_names=None, add_bounding_box_sidesets=False)[source]

Read a mesh file through meshio and convert it to a dualmesh mesh.

Physical groups (Gmsh) or element blocks (Exodus) become subdomain ids, and lower-dimensional cell groups become side sets named after their physical name (or "boundary_<id>" when they are unnamed). boundary_names can rename them, for example {"boundary_1": "inlet"}.

Parameters:

  • filename (str)

  • file_format (str | None)

  • boundary_names (dict | None)

  • add_bounding_box_sidesets (bool)

Return type:

dualmesh._core.Mesh

dualmesh.write_mesh(mesh, filename, file_format=None, **point_data)[source]

Write a mesh (and optional nodal fields) through meshio.

Parameters:
Return type:

None

dualmesh.mesh_from_arrays(points, cells, element_type=None, blocks=None, dimension=None)[source]

Build a mesh from a node array and a connectivity array.

points has shape (num_nodes, 1..3). cells is either an array of shape (num_elements, nodes_per_element) with a single element_type, or a list of (element_type, connectivity) pairs. blocks optionally gives one subdomain id per element.

Parameters:

  • element_type (str | None)

  • dimension (int | None)

Return type:

dualmesh._core.Mesh

dualmesh.graded_coordinates(start, end, num_elements, bias=1.0)[source]

Node coordinates from start to end with a geometric grading.

bias is the ratio between the lengths of successive elements: 1.0 gives a uniform mesh, values larger than one make the elements grow towards end, values smaller than one refine towards end.

Parameters:
Return type:

list[float]

dualmesh.annulus_coordinates(inner_radius, outer_radius, num_elements, bias=1.0)[source]

Radial node coordinates for an annulus (a convenience alias).

Parameters:
Return type:

list[float]

dualmesh.meshing.coordinates_from_spacings(start, spacings)[source]

Node coordinates built by accumulating the given element sizes.

Parameters:
Return type:

list[float]

dualmesh.meshing.mesh_summary(mesh)[source]

Human-readable description of a mesh.

Parameters:

mesh (dualmesh._core.Mesh)

Return type:

str

Results and expressions

class dualmesh.SolveResult(*args, **kwargs)
class dualmesh.ParsedFunction(*args, **kwargs)

Parallel execution

class dualmesh.DistributedProblem(mesh, method='dmcdm', coordinates='cartesian', partitioner='graph', linear_solver='bicgstab', preconditioner='two_level_schwarz', overlap=1, subdomain_solver='ilu', linear_tolerance=1e-10, linear_max_iterations=5000, verbose=False)[source]

A problem split across MPI ranks.

Build it, define the physics on local exactly as for a serial dualmesh.Problem, then call solve() or solve_transient() on this object rather than on the local one.

Parameters:

  • mesh – The whole mesh. Every rank reads it and then keeps only its own part, so the mesh itself is replicated; the degrees of freedom are not.

  • method (str) – As for dualmesh.Problem.

  • coordinates (str) – As for dualmesh.Problem.

  • partitioner (str) – Which partitioner to use; see partition_mesh().

  • linear_solver (str) – "bicgstab" (the default, for any matrix) or "cg" (symmetric positive definite matrices only, about twice as cheap per iteration). Both are distributed: the matrix-vector product multiplies by each rank’s local matrix and adds the results across the partition boundary, and the inner products count every degree of freedom once. The conjugate gradient method needs a symmetric preconditioner, so with it the Schwarz options are applied in their classical, unrestricted form (every subdomain’s whole correction is added) and the coarse level additively; this is weaker than the restricted form, and on the 48 by 48 Poisson problem "cg" takes 64 to 69 iterations on four to eight ranks with any of the three preconditioners.

  • preconditioner (str) –

    "two_level_schwarz", the default, is restricted additive Schwarz on overlapping subdomains with a coarse level. Every rank’s subdomain is its own elements extended by overlap layers of its neighbours’ elements. The subdomain matrix is the global matrix restricted to the subdomain, \(R_\delta A R_\delta^T\), with every row fully assembled (the ranks send each other the entries they hold). Each rank solves its subdomain problem approximately and keeps the correction only on the degrees of freedom it owns (Cai and Sarkis, 1999). The coarse level has one unknown per subdomain and variable (Nicolaides, 1987) and is applied before the subdomain solves (the operator called A-DEF1 by Tang, Nabben, Vuik and Erlangga, 2009), so that information crosses the whole mesh in one application.

    "additive_schwarz" is the same without the coarse level. "jacobi" divides by the diagonal of the global matrix.

    On the Poisson problem of a 48 by 48 mesh, BiCGSTAB to a relative tolerance of \(10^{-10}\) takes the following numbers of iterations (measured with this library):

    ranks

    1

    2

    4

    8

    16

    two_level_schwarz

    24

    25

    23

    25

    24

    additive_schwarz

    24

    26

    29

    33

    33

    jacobi

    50

    49

    49

    50

    49

    On a 128 by 128 mesh the two-level counts are 66, 60 and 64 on one, four and sixteen ranks, the one-level counts 66, 76 and 86, and Jacobi’s 137, 131 and 132. The iteration count is reported as linear_iterations on the result of solve(), and it is the number to watch: a preconditioner whose count grows in proportion to the number of ranks cancels the benefit of the extra ranks.

  • overlap (int) – Layers of elements by which each Schwarz subdomain reaches into its neighbours (default 1). Zero gives non-overlapping subdomains. More overlap lowers the iteration count and raises the cost of each subdomain solve and of the setup exchange.

  • subdomain_solver (str) – "ilu" (the default: incomplete LU without fill, cheap to build and apply) or "lu" (an exact sparse LU of each subdomain matrix, which roughly halves the iteration count in two dimensions and is expensive for large three-dimensional subdomains).

  • linear_tolerance (float) – As for the serial solver.

  • linear_max_iterations (int) – As for the serial solver.

  • verbose (bool) – As for the serial solver.

local

The rank-local problem; define the physics on it.

property rank: int

This process’s rank.

property num_ranks: int

The number of processes the problem is split across.

property num_owned_dofs: int

Degrees of freedom this rank owns.

property num_global_dofs: int

Degrees of freedom of the whole problem.

solve(**options)[source]

Solve the steady problem. Takes the same options as dualmesh.Problem.solve(), except that the linear solver is configured in the constructor.

solve_transient(end_time, dt, start_time=0.0, theta=1.0, output_interval=0, output_file_base='', time_stepper='fixed', dt_min=0.0, dt_max=0.0, growth_factor=2.0, cutback_factor=0.5, error_tolerance=0.001, optimal_iterations=4, iteration_window=2, max_rejected_steps=10, **options)[source]

Advance the problem in time, with the same arguments and the same time steppers as dualmesh.Problem.solve_transient().

The adaptive controllers are driven by quantities that are already global reductions – the residual norm and the norm of the difference between the coarse and the fine step – so every rank reaches the same decision and they step in lockstep without any extra communication. With output_file_base set, each output step writes one .vtu per rank and one .pvtu index.

Parameters:
gathered_values(variable)[source]

Values of a variable at every node of the whole mesh, assembled on every rank.

This allocates one vector of the global size per rank, so it is meant for testing and for small problems; large runs should write the result with write_vtu() instead.

Parameters:

variable (str)

write_vtu(base, cell_properties=())[source]

Write one .vtu file per rank plus a .pvtu index.

base is the file name without an extension. ParaView and VisIt open the .pvtu file as a single data set.

Parameters:
Return type:

None

summary()[source]

One line describing the partition, for logs.

Return type:

str

dualmesh.partition_mesh(mesh, num_parts, method='recursive_coordinate_bisection')[source]

Split a mesh into num_parts groups of elements.

Parameters:

  • mesh – The mesh to split.

  • num_parts (int) – How many parts to make. It must not exceed the number of elements.

  • method (str) – "recursive_coordinate_bisection" repeatedly halves the set of element centroids along its longest axis. It is fast and deterministic and needs no connectivity, but it cuts more faces than necessary on an unstructured mesh. "graph" grows each part outward from a seed element through the face connectivity, so the parts follow the mesh topology. "metis" calls METIS on the dual graph of the mesh and normally gives the smallest cut; it falls back to "graph" when the extension was built without METIS.

Returns:

element_part (the part of every element), node_owner (the part that owns every node, which is the smallest part index among those that touch it), edge_cut (the number of faces whose two elements are in different parts, the quantity a partitioner tries to minimise), and largest_part and smallest_part (element counts, a measure of load balance).

Return type:

dict

dualmesh.have_mpi()[source]

Whether the extension was built with MPI support.

Return type:

bool

dualmesh.have_metis()[source]

Whether the extension was built against METIS.

METIS usually cuts fewer faces than the built-in partitioners on unstructured meshes. It is the partitioner libMesh, and therefore MOOSE, uses.

Return type:

bool

dualmesh.is_root()[source]

Whether this is rank 0, the process that should do the printing.

Return type:

bool

dualmesh.num_ranks()[source]

The number of processes; 1 in a serial run.

Return type:

int

Adaptive refinement

dualmesh.solve_with_adaptive_refinement(build_problem, mesh, variable, num_cycles=3, marker=<function mark_by_error_fraction>, max_elements=None, callback=None)[source]

Solve, estimate, mark and refine, num_cycles times.

build_problem is called with a mesh and must return a solved-ready Problem defined on it: the whole problem is rebuilt on every mesh rather than transferred, which keeps the boundary conditions and the material state exactly as the user wrote them. marker turns the vector of indicators into a boolean array of elements to refine. max_elements stops the loop early once the mesh has grown past that size, and callback(cycle, problem, indicators) is called after every solve, which is where a calculation writes output or records a convergence history.

Returns (problem, mesh) for the last, finest mesh.

Parameters:
dualmesh.refine_marked(mesh, marked)[source]

Refine the marked elements, keeping the mesh conforming.

marked holds one flag per element. More elements than those marked are bisected, because the neighbours of a bisected edge have to be bisected too. Returns (refined_mesh, parents), where parents[i] is the index in mesh of the element that element i of the refined mesh came from; a field stored per element is carried over with it.

Parameters:

marked (Sequence[bool])

dualmesh.mark_by_fraction(indicators, fraction=0.3)[source]

Mark the given fraction of the elements with the largest indicators.

This is the simplest marking rule: with fraction=0.3 the worst thirty per cent of the elements are refined at every cycle, so the mesh grows at a predictable rate whatever the shape of the error distribution.

Parameters:

fraction (float)

Return type:

numpy.ndarray

dualmesh.mark_by_error_fraction(indicators, fraction=0.5)[source]

Mark the smallest set of elements that carries the given share of the total squared error.

This is bulk, or Dörfler, marking: with fraction=0.5 just enough elements are refined to account for half of the total error. It refines few elements when the error is concentrated in a few of them and many when the error is spread out, which is the behaviour that proofs of convergence for adaptive methods rely on.

W. Dörfler, “A convergent adaptive algorithm for Poisson’s equation”, SIAM Journal on Numerical Analysis 33 (1996) 1106-1124.

Parameters:

fraction (float)

Return type:

numpy.ndarray

dualmesh.mark_by_threshold(indicators, threshold)[source]

Mark every element whose indicator exceeds an absolute threshold.

Parameters:

threshold (float)

Return type:

numpy.ndarray

Manufactured solutions

Verification by the method of manufactured solutions.

Background

A code is verified when it is shown to solve its equations correctly, which is a different question from whether the equations describe nature. The standard tool is the method of manufactured solutions (Roache, 2002): choose any smooth function as the exact solution, substitute it into the governing equations, and whatever is left over is the source that function needs. Adding that source to the problem and imposing the chosen function on the boundary gives a problem whose exact solution is known, however nonlinear or coupled the equations are. Solving it on a sequence of refined meshes then measures the order at which the error falls, and a discretisation that is correctly implemented must show its theoretical order. An error in the implementation almost always shows up as a lower order, often as no convergence at all.

The weak point of the method is keeping the two descriptions of the equations in step: the one the code solves, and the one the source is derived from. A hand-derived source for a nonlinear, axisymmetric, coupled problem is itself an invitation to error. This module removes that step. Each Term knows both which dualmesh objects it adds and, in SymPy, which flux and source those objects contribute, so the source is derived from exactly the equations the code assembles, by the same canonical form

\[\mathcal{R}(u) = -\nabla \cdot \mathbf{F}(u, \nabla u, \mathbf{x}, t) + S(u, \nabla u, \mathbf{x}, t) = 0 ,\]

with the divergence taken in the problem’s coordinate system.

Example

from dualmesh import mms

study = mms.ManufacturedSolution(
    fields={"u": "sin(pi*x)*cos(pi*y)"},
    terms=[mms.Diffusion("u", diffusivity="1 + x*y", polynomial=(1, 0.5)),
           mms.Advection("u", velocity=(1, 0.5, 0))],
    dimension=2,
)
result = study.convergence_study(
    lambda n: dm.generate_rectangle_mesh(x_min=0, x_max=1, y_min=0, y_max=1,
                                         num_x_elements=n, num_y_elements=n),
    levels=[4, 8, 16, 32], method="dmcdm")
print(result.table())
result.rates("u", "l2")      # tends to 2 for linear elements

SymPy is needed only for this module, which imports it lazily.

References: P. J. Roache, “Code verification by the method of manufactured solutions”, Journal of Fluids Engineering 124 (2002) 4-10; K. Salari and P. Knupp, Code Verification by the Method of Manufactured Solutions, Sandia Report SAND2000-1444, 2000.

class dualmesh.mms.Advection(variable, velocity=(1.0, 0.0, 0.0), form='non_conservative')[source]

Advection kernel with a constant velocity. The non-conservative form contributes the source \(\mathbf{v} \cdot \nabla u\), the conservative form the flux \(-\mathbf{v} u\).

Parameters:
contributions(fields, coordinates, dimension)[source]

Return {variable: (flux, source)} for the exact fields, with the flux a list of three SymPy expressions.

class dualmesh.mms.ConvergenceResult(method, sizes=<factory>, num_dofs=<factory>, errors=<factory>)[source]

Errors on a sequence of meshes, and the observed orders.

sizes holds the characteristic element size \(h\) of every mesh, taken as \((|\Omega| / N_e)^{1/d}\), which reduces to the element size for a uniform mesh and is well defined for any mesh.

Parameters:
rates(variable, norm='l2')[source]

Observed order between consecutive meshes, \(\log(e_{k}/e_{k+1}) / \log(h_{k}/h_{k+1})\).

Parameters:
Return type:

list

order(variable, norm='l2', last=2)[source]

Least-squares slope of \(\log e\) against \(\log h\) over the last finest meshes, which is the asymptotic order to report.

Parameters:
Return type:

float

class dualmesh.mms.Diffusion(variable, diffusivity=1.0, polynomial=(1.0,))[source]

Diffusion kernel: \(\mathbf{F} = k(\mathbf{x}, t)\, p(u)\, \nabla u\) with \(p(u) = c_0 + c_1 u + c_2 u^2 + \dots\).

A non-constant polynomial makes the problem nonlinear.

Parameters:

  • variable (str)

  • polynomial (Sequence[float])

contributions(fields, coordinates, dimension)[source]

Return {variable: (flux, source)} for the exact fields, with the flux a list of three SymPy expressions.

class dualmesh.mms.LinearElasticity(displacements, youngs_modulus, poissons_ratio, formulation='plane_strain')[source]

LinearElasticStress material and one StressDivergence kernel per displacement component. The flux of component \(i\) is row \(i\) of the stress; in the axisymmetric formulation the radial equation also carries the hoop source \(\sigma_{\theta\theta} / r\).

Parameters:

  • displacements (Sequence[str])

  • youngs_modulus (float)

  • poissons_ratio (float)

  • formulation (str)

contributions(fields, coordinates, dimension)[source]

Return {variable: (flux, source)} for the exact fields, with the flux a list of three SymPy expressions.

class dualmesh.mms.ManufacturedSolution(fields, terms, dimension, coordinates='cartesian')[source]

A problem whose exact solution is chosen, not computed.

fields maps each variable name to its exact solution, as text or as a SymPy expression in x, y, z and t. terms are the pieces of the governing equations. coordinates is the problem’s coordinate system, which decides the form of the divergence:

  • Cartesian: \(\nabla \cdot \mathbf{F} = \sum_d \partial F_d / \partial x_d\);

  • axisymmetric, with \(x = r\) and \(y = z\): \(r^{-1} \partial (r F_r) / \partial r + \partial F_z / \partial z\);

  • spherical, with \(x = r\): \(r^{-2} \partial (r^2 F_r) / \partial r\).

Parameters:

  • fields (dict)

  • terms (Sequence[Term])

  • dimension (int)

  • coordinates (str)

forcing()[source]

The source each equation needs for the chosen fields to satisfy it: \(f = -\nabla \cdot \mathbf{F}(u) + S(u)\), which a BodyForce of intensity \(f\) (whose own source is \(-f\)) cancels.

Return type:

dict

build(mesh, method='dmcdm', boundary=None, start_time=0.0, **problem_options)[source]

A dualmesh problem for this manufactured solution on mesh: the variables, the terms, one BodyForce per equation carrying the manufactured source, and the exact solution imposed on boundary (every side set when omitted). The unknowns start at the exact solution at start_time, which is the initial condition of a transient study and a good first iterate for a nonlinear one.

Parameters:
Return type:

Problem

errors(problem)[source]

{(variable, "l2"): value, (variable, "h1"): value} for a solved problem, at the problem’s current time.

Parameters:

problem (Problem)

Return type:

dict

convergence_study(mesh_factory, levels, method='dmcdm', transient=None, solve_options=None, **problem_options)[source]

Solve on mesh_factory(n) for every n in levels and record the errors.

transient, when given, holds the arguments of solve_transient(); its dt may be a callable of the element size \(h\), so that the time step is refined together with the mesh, which is what a study of the combined space-time order needs.

Parameters:
Return type:

ConvergenceResult

class dualmesh.mms.Reaction(variable, coefficient=1.0, exponent=1.0)[source]

Reaction kernel: the source \(c\, u^p\).

Parameters:
contributions(fields, coordinates, dimension)[source]

Return {variable: (flux, source)} for the exact fields, with the flux a list of three SymPy expressions.

class dualmesh.mms.Term[source]

One piece of the governing equations.

A term does two things, and the point of the class is that the two cannot drift apart: add_to() adds the dualmesh objects, and contributions() returns, in SymPy, the flux and the source those objects put into the canonical form.

contributions(fields, coordinates, dimension)[source]

Return {variable: (flux, source)} for the exact fields, with the flux a list of three SymPy expressions.

Parameters:
Return type:

dict

class dualmesh.mms.TimeDerivative(variable, coefficient=1.0)[source]

TimeDerivative kernel: the source \(c\, \partial u / \partial t\) of the continuous equation.

Parameters:

variable (str)

contributions(fields, coordinates, dimension)[source]

Return {variable: (flux, source)} for the exact fields, with the flux a list of three SymPy expressions.

Physics helpers

Ready-made physics: one call adds all the kernels of a model.

The models of the solid mechanics and fluids modules need one kernel per equation (and, for shear-deformable plates, a second kernel for the transverse shear terms so that they can be integrated with a reduced rule). These helpers add them consistently, in the spirit of MOOSE’s Physics syntax:

import dualmesh as dm

problem = dm.Problem(mesh)
dm.physics.add_plane_elasticity(
    problem, youngs_modulus=30e6, poissons_ratio=0.25, thickness=0.036)
dualmesh.physics.add_plane_elasticity(problem, displacements=('displacement_x', 'displacement_y'), youngs_modulus=1.0, poissons_ratio=0.0, formulation='plane_stress', thickness=1.0, body_force=None, name='elasticity', **material_parameters)[source]

Add the two (or three) equilibrium equations of linear elasticity.

Creates the displacement variables if they do not exist yet, a LinearElasticStress material, and one StressDivergence kernel per component. formulation is "plane_stress", "plane_strain", "axisymmetric", or "three_dimensional".

Parameters:
dualmesh.physics.add_incompressible_flow(problem, velocities=('velocity_x', 'velocity_y'), dynamic_viscosity=1.0, density=0.0, penalty_parameter=100000000.0, body_force=None, name='flow')[source]

Add the penalty momentum equations of an incompressible flow.

With density = 0 the equations are the Stokes equations; a positive density adds the convective term of the Navier-Stokes equations. A PenaltyPressure material provides the recovered pressure as the material property "pressure".

Parameters:
dualmesh.physics.add_boussinesq_buoyancy(problem, temperature, velocities, gravity, density=1.0, thermal_expansion=1.0, reference_temperature=0.0, scale_with_load=False, name='buoyancy')[source]

Add the Boussinesq buoyancy force to the momentum equations.

One BoussinesqBuoyancy kernel is added for every velocity component along which gravity acts, giving the body force \(\mathbf{f} = -\rho_0 \beta (T - T_0) \mathbf{g}\). With scale_with_load=True load stepping ramps the buoyancy, which is how a high Rayleigh number is reached from rest. Together with a HeatConvection kernel on the temperature this makes a natural convection problem; see examples/natural_convection.py.

Parameters:
dualmesh.physics.add_beam(problem, model='BeamEulerBernoulliMixed', axial_displacement='axial_displacement', transverse_displacement='deflection', third_variable=None, name='beam', **parameters)[source]

Add a beam model (three kernels, one per variable).

model is "BeamEulerBernoulliMixed", "BeamTimoshenkoMixed", or "BeamTimoshenkoDisplacement". For the displacement Timoshenko model the shear term of the rotation equation is integrated with a reduced rule, which is what prevents shear locking.

Parameters:

  • model (str)

  • axial_displacement (str)

  • transverse_displacement (str)

  • third_variable (str | None)

  • name (str)

dualmesh.physics.add_circular_plate(problem, radial_displacement='radial_displacement', transverse_displacement='deflection', rotation='rotation', bending_moment='bending_moment', theory='first_order', name='plate', **parameters)[source]

Add an axisymmetric circular plate model on a radial mesh.

The problem must use coordinates="axisymmetric", so that every integral carries the factor \(2 \pi r\).

With theory="first_order" (the default) this adds the first-order shear deformation, or Mindlin, model in terms of the radial displacement \(u\), the deflection \(w\), and the rotation \(\phi_r\). Two kernels are added per variable: one for the bending and membrane terms, and one for the transverse shear force, which is integrated at the centre of the element so that thin plates do not lock.

With theory="classical" this adds the mixed classical, or Kirchhoff, model in terms of \(u\), \(w\), and the radial bending moment \(M_{rr}\). The classical theory gives a fourth-order equation in \(w\), which the dual mesh control domain method cannot discretize, so the bending moment is carried as a third unknown and the system becomes three second-order equations. There is no shear term to under-integrate, so one kernel per variable is enough. The natural boundary quantity of the moment equation is the slope \(dw/dr\), which means a clamped edge needs no condition at all on that equation, while a simply supported edge prescribes \(M_{rr} = 0\).

Parameters:

  • radial_displacement (str)

  • transverse_displacement (str)

  • rotation (str)

  • bending_moment (str)

  • theory (str)

  • name (str)

dualmesh.physics.add_plate(problem, in_plane_displacements=('displacement_x', 'displacement_y'), transverse_displacement='deflection', rotations=('rotation_x', 'rotation_y'), name='plate', **parameters)[source]

Add the first-order (Mindlin) rectangular plate model (five variables).

Two kernels are added per variable: one for the bending and membrane terms and one for the transverse shear forces, the latter with reduced integration, which removes shear locking in thin plates.

Parameters:

Objects written in Python

Base classes for kernels, boundary conditions, and materials written in Python.

Subclass these to add physics without touching C++. Inside the callbacks the solution is delivered as dualmesh.ADReal numbers, so the Jacobian used by Newton’s method remains exact:

import dualmesh as dm

class NonlinearBar(dm.PythonKernel):
    "Axial bar with a = EA (1 + 1.5 u' + 0.5 u'^2)."

    def setup(self, problem):
        self.axial_stiffness = self.parameters["axial_stiffness"]

    def has_flux(self):
        return True

    def compute_flux(self, ctx):
        strain = ctx.gradient(self.variable)[0]
        a = self.axial_stiffness * (1.0 + 1.5 * strain + 0.5 * strain * strain)
        return [a * strain]
class dualmesh.objects.PythonKernel(*args, **kwargs)[source]

Bases: Kernel, _PythonObjectMixin

A kernel implemented in Python (flux and/or source).

Parameters:

  • variable (str)

  • name (str | None)

  • quadrature (str)

  • reduced_integration (bool)

  • scale_with_load (bool)

class dualmesh.objects.PythonBoundaryCondition(*args, **kwargs)[source]

Bases: IntegratedBC, _PythonObjectMixin

An integrated boundary condition implemented in Python.

compute_boundary_flux returns the outward normal flux q = n . F.

Parameters:

  • variable (str)

  • name (str | None)

  • quadrature (str)

  • reduced_integration (bool)

  • scale_with_load (bool)

class dualmesh.objects.PythonNodalBoundaryCondition(*args, **kwargs)[source]

Bases: NodalBC, _PythonObjectMixin

An essential (Dirichlet) boundary condition implemented in Python.

Parameters:

  • variable (str)

  • name (str | None)

  • scale_with_load (bool)

class dualmesh.objects.PythonMaterial(*args, **kwargs)[source]

Bases: Material, _PythonObjectMixin

A material implemented in Python.

Declare properties in declare_properties and fill them in compute_properties:

class TemperatureDependentConductivity(dm.PythonMaterial):
    def declare_properties(self, registry):
        self.conductivity_id = registry.declare("thermal_conductivity")

    def setup(self, problem):
        self.temperature = problem.variable_index("temperature")

    def compute_properties(self, ctx):
        T = ctx.coefficient_value(self.temperature)
        ctx.set_property(self.conductivity_id, 20.0 + 0.2 * T)
Parameters:

name (str | None)

setup(problem)[source]

Resolve names (variables, functions, properties) before solving.

Functionally graded sections

Stiffness coefficients of functionally graded beams and plates.

A functionally graded material (FGM) varies continuously through the thickness. Following Reddy, the modulus follows the power law

\[E(z) = (E_1 - E_2) \left( \frac{1}{2} + \frac{z}{h} \right)^{n} + E_2 , \qquad -\frac{h}{2} \le z \le \frac{h}{2},\]

so that \(E(h/2) = E_1\) (top) and \(E(-h/2) = E_2\) (bottom); the power-law index \(n = 0\) gives a homogeneous beam with \(E = E_1\).

The stress resultants of a beam of width \(b\) use

\[(A_{xx}, B_{xx}, D_{xx}) = b \int_{-h/2}^{h/2} E(z) (1, z, z^2) \, dz , \qquad S_{xz} = \frac{K_s}{2(1+\nu)} b \int_{-h/2}^{h/2} E(z) \, dz ,\]

with \(K_s = 5/6\) the shear correction factor of a rectangular section. For plates every coefficient is divided by \(1 - \nu^2\) (plane-stress reduced stiffness) and \(b = 1\) is taken per unit width.

class dualmesh.fgm.BeamStiffness(extensional, coupling, bending, shear)[source]

Extensional, coupling, bending, and shear stiffnesses of a beam.

Parameters:
property reduced_bending: float

The reduced bending stiffness \(D^{*} = D_{xx} A_{xx} - B_{xx}^2\).

property effective_extensional: float

Reddy’s \(\bar{A}_{xx} = D^{*} / D_{xx}\) (mixed formulations).

property effective_coupling: float

Reddy’s \(\bar{B}_{xx} = B_{xx} / D_{xx}\) (mixed formulations).

dualmesh.fgm.beam_stiffness(modulus_top, modulus_bottom, power_law_index, height, width=1.0, poisson_ratio=0.3, shear_correction_factor=0.8333333333333334, plate=False)[source]

Closed-form stiffnesses of a power-law functionally graded section.

modulus_top is \(E_1\) (at \(z = +h/2\)) and modulus_bottom is \(E_2\) (at \(z = -h/2\)). With plate=True the coefficients are divided by \(1 - \nu^2\), as required by the plate theories.

Parameters:
Return type:

BeamStiffness

dualmesh.fgm.modulus(z, modulus_top, modulus_bottom, power_law_index, height)[source]

The through-thickness modulus \(E(z)\) of the power-law profile.

Parameters:
Return type:

float

dualmesh.fgm.stiffness_by_quadrature(modulus_top, modulus_bottom, power_law_index, height, width=1.0, num_points=400)[source]

(A, B, D) obtained by numerical integration (used to check the formulas).

Parameters:
Return type:

tuple[float, float, float]

Expressions

Expressions given as text, turned into functions of (x, y, z, t).

This is the counterpart of MOOSE’s ParsedFunction. A coefficient, a source or a boundary value can be written as an expression:

import dualmesh as dm

top = dm.parsed_function("500*(1 - 10*x^2)")
problem.add_boundary_condition(
    "DirichletBC", variable="temperature", boundary="top", value=top)

or, more briefly, passed as the text itself, which is compiled the same way:

problem.add_kernel("BodyForce", variable="u", value="sin(pi*x)*exp(-t)")

The expression is compiled once, to a short program evaluated in C++. That matters for speed and for threading: a Python callable would be called back at every quadrature point of every element on every iteration, and because calling into Python needs the interpreter lock it would also force the assembly onto a single thread. A compiled expression does neither.

The grammar is the usual one. + - * / have their usual precedence, ^ and ** both mean exponentiation and bind tighter than a unary minus, and comparisons < > <= >= == != give 1 or 0. The names in scope are x, y, z, t, the constants pi and e, and the functions sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, exp, log, log10, sqrt, abs, floor, ceil, erf, sign, atan2, pow, hypot, min, max and if(condition, a, b). SymPy’s printed form (E, Abs, **) is accepted unchanged, so a manufactured source derived symbolically can be passed straight in.

dualmesh.expressions.parsed_function(expression)[source]

Compile an expression in x, y, z and t.

The result is a dualmesh.ParsedFunction, which every parameter that takes a function accepts, and which can also be called from Python as f(x, y, z, t) to check it. A syntax error is reported with its position in the text.

Parameters:

expression (str)

Return type:

dualmesh._core.ParsedFunction

Post-processing

Small helpers for extracting and comparing results.

dualmesh.postprocess.sample_line(problem, variable, start, end, num_points=21)[source]

Sample a variable along a straight line; returns (arc_length, values).

Parameters:

  • variable (str)

  • num_points (int)

dualmesh.postprocess.values_on_line(problem, variable, coordinates, axis=0, other=0.0)[source]

Sample a variable at given coordinates along one axis.

Parameters:
dualmesh.postprocess.relative_error(computed, reference)[source]

Element-wise relative error, using the reference magnitude as the scale.

Return type:

numpy.ndarray

dualmesh.postprocess.convergence_rates(mesh_sizes, errors)[source]

Observed convergence rates between successive refinements.

Parameters:
Return type:

numpy.ndarray

dualmesh.postprocess.comparison_table(labels, computed, reference, title='')[source]

A fixed-width table of computed values against reference values.

Parameters:
Return type:

str

Automatic differentiation

Elementary functions that work for both floats and AD numbers.

Kernels written in Python receive dualmesh.ADReal values, whose derivatives with respect to the local degrees of freedom are carried along so that the Jacobian used by Newton’s method stays exact. Use these functions instead of math inside such kernels:

import dualmesh as dm

class ArrheniusReaction(dm.PythonKernel):
    def compute_source(self, ctx):
        temperature = ctx.value(self.temperature_index)
        return self.pre_exponential * dm.exp(-self.activation / temperature)
dualmesh.ad.sqrt(value)

sqrt(x) for floats and AD numbers.

dualmesh.ad.exp(value)

exp(x) for floats and AD numbers.

dualmesh.ad.log(value)

log(x) for floats and AD numbers.

dualmesh.ad.sin(value)

sin(x) for floats and AD numbers.

dualmesh.ad.cos(value)

cos(x) for floats and AD numbers.

dualmesh.ad.tanh(value)

tanh(x) for floats and AD numbers.

dualmesh.ad.abs(value)[source]

Absolute value for floats and AD numbers.

dualmesh.ad.pow(base, exponent)[source]

Power for floats and AD numbers.

Command line

Command-line driver: run a problem described in an input file.

Input files are YAML and mirror the block structure of the Python API (and of MOOSE input files):

mesh:
  type: rectangle
  x_min: 0.0
  x_max: 0.1
  y_min: 0.0
  y_max: 0.05
  num_x_elements: 10
  num_y_elements: 5

problem:
  method: dmcdm            # or fem, hfvm, zfvm
  coordinates: cartesian   # or axisymmetric, spherical
  threads: 4               # optional; the default uses every core

variables:
  temperature: {initial_condition: 0.0}

functions:
  ambient: "40 + 10*x"     # an expression in x, y, z, t

kernels:
  conduction:
    type: HeatConduction
    variable: temperature
    thermal_conductivity: 20.0
  heating:
    type: HeatSource
    variable: temperature
    heat_source: 1.0e6

boundary_conditions:
  left:  {type: DirichletBC, variable: temperature, boundary: left, value: 40}
  right: {type: DirichletBC, variable: temperature, boundary: right, value: 10}
  top:
    type: ConvectiveHeatFluxBC
    variable: temperature
    boundary: top
    heat_transfer_coefficient: 75.0

executioner:
  type: steady             # or transient
  nonlinear_solver: newton

outputs:
  vtu: bus_bar.vtu
  csv: bus_bar.csv
  reactions: [[temperature, left]]

Usage:

dualmesh run input.yaml
mpirun -n 4 dualmesh run input.yaml     # distributed, when built with MPI
dualmesh list --category Kernel
dualmesh describe HeatConduction

Under mpirun with more than one process the problem is solved by DistributedProblem, configured by an optional parallel block (partitioner, linear_solver, preconditioner, overlap, subdomain_solver, linear_tolerance, linear_max_iterations); the same input file runs unchanged on one process.

dualmesh.cli.build_mesh(block)[source]

Build a mesh from the mesh block of an input file.

Parameters:

block (dict[str, Any])

dualmesh.cli.build_problem(document)[source]

Build the problem described by a parsed input file.

Returns the Problem, or, when the program runs on more than one MPI process, the DistributedProblem whose local problem carries the objects.

Parameters:

document (dict[str, Any])

dualmesh.cli.run(document, verbose=False)[source]

Build and solve the problem described by an input file.

Parameters: