Poly-memb

poly-memb is a Python library implementing a 2D simulator of dynamics of membranes subject to electric fields. The module ‘mesh_manager’ provides functions to generate polytopal meshes by cutting and agglomerating background meshes intersected by moving interfaces. The module ‘solvers’ provides functions to implement polygonal numerical schemes. The simulator is based on a Hybrid High Order scheme [1] to solve flow variables and a Discrete De Rham [2] scheme to solve for electric variables.

Gallery

Example of supported mesh, before and after agglomeration. The polygonal interface cuts originally simplicial elements to originate conformal polygonal subelements. The cut procedure is reiterated after movement of the membrane.

Electrostatic potential and electric field resulting from the DDR method. The conformal setting allows for seamless treatment of interface conditions.

A HHO solver is used to determine the membrane velocity. Interface edges are embedded as mesh edges.

Modules

Module: mesh_manager

polymemb.mesh_manager.R2_norm(x)

Calculates Euclidean norm in R2

Parameters:

x (np.ndarray) – R2 vector whose to calculate norm of

Returns:

float

polymemb.mesh_manager.adjacent_elements(elem2node, cut_iel, intsec_points)

List of adjacent elements (indexes)

Parameters:

• elem2node (list(list(in))) – elem2node connectivity

• cut_iel (int) – element index

Returns:

adjacent elements

Return type:

list(int)

polymemb.mesh_manager.agglomerate(mesh, ref_mesh, cryt_size, cryt_skewness, verbose=False)

Agglomerate bad quality elements :type mesh: mema.mesh2D :param mesh: mesh :type mesh: mema.mesh2D :type ref_mesh: mema.mesh2D :param ref_mesh: reference uncut mesh :type ref_mesh: mema.mesh2D :type cryt_size: float :param cryt_size: minimal accepted relative size for elements :type cryt_size: float :type cryt_skewness: float :param cryt_skewness: maximal accepted skewness for elements :type cryt_skewness: float :type verbose: boolean :param verbose: display partial outputs :type verbose: boolean

polymemb.mesh_manager.agglomerate_2_elems(elem_agglomerating, elem_to_agglomerate, intface_data_agglomerating, intface_data_to_agglomerate, verbose=False)

Agglomerate 2 elements provided as lists of points ALERT: currently, agglomerations are done under hp that elements share a connected chain of points

Args: elem_agglomerating (list(int)): agglomerating element as list of points elem_to_agglomerate (list(int)): element to agglomerate as list of points intface_data_agglomerating (list(list(int))): interface data (formatted as detailed in get_intface_data) for agglomerating element intface_data_to_agglomerate(list(list(int))): interface data (formatted as detailed in get_intface_data) for element to agglomerate verbose (boolean): display in-function outputs

polymemb.mesh_manager.barycenter(polygon)

Compute baricenter of a polygon

Args: polygon (list(np.array))

polymemb.mesh_manager.bm(mesh, intface, execute_fill_side_mask=True, verbose=False)

Breaks mesh along interface Remark: treat also case of element with multiple cuts

Parameters:

• mesh (mesh2D) – mesh to break

• intface (disk_interface) – interface along which to cut

Returns:

new broken mesh, with list of gamma_couples, mesh2intface

and side mask

Return type:

mesh2D

polymemb.mesh_manager.break_mesh(mesh, intface, execute_fill_side_mask=True, verbose=False)

Breaks mesh along interface Remark: treat also case of element with multiple cuts

Parameters:

• mesh (mesh2D) – mesh to break

• intface (disk_interface) – interface along which to cut

Returns:

new broken mesh, with list of gamma_couples, mesh2intface

and side mask

Return type:

mesh2D

polymemb.mesh_manager.calc_intersection(x1, x2, y1, y2, rho)

Check if segment on a plane intersects with ball of radius rho

Parameters:

• x1 (float) – x coord of point1

• x2 (float) – x coord of point2

• y1 (float) – y coord of point1

• y2 (float) – y coord of point2

• rho (float) – disk radius

Returns:

segment is intersected numpy.ndarray: intersection position

Return type:

boolean

polymemb.mesh_manager.calc_intersection_segments(A1, B1, A2, B2)

polymemb.mesh_manager.calc_membrane_energy(intface)

Calculates total energy of the membrane

polymemb.mesh_manager.calc_membrane_power(intface)

Calculates power acting on the membrane

polymemb.mesh_manager.calc_membrane_volume(intface)

Calculates volume of surface

Parameters:

u (np.array) – P^1(Gamma_h)

polymemb.mesh_manager.calc_time_step(intface, eta)

Calculates time step such that energy increase is controlled :type intface: mema.disk_intface :param intface: interface :type intface: mema.disk_intface :type eta: float :param eta: user-dependent parameter :type eta: float

polymemb.mesh_manager.check_and_add_new_points(new_coords, new_points)

Add a bunch of new points to a list of points if not already there and provide indexes on the new entries

Parameters:

• new_coords (list(np.array)) – list of points

• new_points (list(np.arra)) – points to add

polymemb.mesh_manager.check_if_in_triangle(p1, p2, p3, q)

Check if point q is inside the triangle formed by points p1, p2, and p3.

Parameters: p1, p2, p3: Tuples representing the vertices of the triangle (x, y). q: Tuple representing the point to check (x, y).

Returns: bool: True if q is inside the triangle, False otherwise.

polymemb.mesh_manager.check_if_present(mesh, node_coords)

polymemb.mesh_manager.compute_aspect_ratio(polygon_coords)

Compute the robust aspect ratio of a polygon using PCA (bounding box in the system of princiapl axes)

polymemb.mesh_manager.cut_elem(coords, elem, intface, ied, initial_enter_point, idx_initial_enter_point, verbose)

Function to cut an element, and enrich a list of mesh points

Parameters:

• coords (list(np.array)) – list of points

• elem (list(int)) – list of node idxs of the element

• intface (disk_interface) – interface

• ied (int) – entering intface edge

• initial_enter_point (np.array) – coords of initial point of break mesh

• idx_initial_enter_points (int) – its index

Returns:

enriched list points list (int): list of intface edges of the cut list(int): first child element list(nt): second child element int: last edge of the cut (and first of the next)

Return type:

list (np.array)

polymemb.mesh_manager.default_cryterion(iel, elem_coords)

Check if cryterion to refine is respected XXX dummy cryterion

Parameters:

• iel (int) – element index

• elem_coords – element node coordinates

Returns:

bool

class polymemb.mesh_manager.disk_interface(edges, coords, k_b, k_str, initial=True, init_edge_length=None, velocity=None)

Bases: object

Interface class; implemented as chain of edges

edges

edge to node connectivity

Type:

list(list(int))

coords

node coordinates

Type:

list(np.array)

velocity

edge velocity

Type:

list(np.array)

k_b

bending modulus

Type:

float

k_str

stretching modulus

Type:

float

init_edge_length

reference edge length

Type:

list(float)

velocity

edge velocity (0 if not assigned)

Type:

list(np.array)

advect(ref_vol, delta_t, vol_corr=True)

Advect interface using velocity

Parameters:

• ref_vol (float) – reference volume

• delta_t (float) – time step

Returns:

moved interface

Return type:

disk_intface

apply_LB(data)

Applies discrete Laplace-Beltrami Operator

Parameters:

data (list) – node-valued field

Returns

list: node_valued curve Laplacian (LB) of data

calc_F_bending()

Calculates nodal force $F_{bending}$;

Returns:

nodal force (vector)

Return type:

list(np.array)

calc_F_stretching(drop_model=False, gamma=1)

Calculates nodal force $F_{stretching}$; Hooke’s elastic law

Returns:

nodal force (vector)

Return type:

list(np.array)

calc_curvature()

Calculate node curvature

Returns:

node curvature

Return type:

list(float)

calc_edge_length()

Calculates length of the edges

Returns:

list of edge lengths

Return type:

list(float)

calc_nodal_forces()

Calculates nodal forces; $F = F_{bending} + F_{stretching}$

Returns:

edge_wise constant surface tension (vector)

Return type:

list(np.array)

calc_normal()

Calculate edge normal

calc_t_gamma()

Calculates surface tension $t_Gamma$ by calculating nodal forces and transferring onto edges with correct scaling

Returns:

edge_wise constant surface tension (vector)

Return type:

list(np.array)

transfer_edge2node(data)

Tranfers an interface property from edges to nodes or viceversa applying average

Parameters:

data (list) – a data defined as a list over edges or nodes

Returns:

ransferred data

Return type:

list

polymemb.mesh_manager.display_element(mesh, elem, ax, color, width=0.0001, display_nodes=True)

Display single element

polymemb.mesh_manager.extract_first_vertex_coordinates(stl_file_path)

Extract list of coordinated from file

Parameters:

stl_file_path (string) – path to file

Returns:

list

polymemb.mesh_manager.fill_side_mask(mesh, no_initial_points, only_initialize=False)

Complete side mask by extrapolating value at interface up to the border (uses dofs associated to nodes as in ddrin)

Parameters:

• mesh (mesh2D) – mesh

• no_initial_points (int) – number of points before break_mesh

polymemb.mesh_manager.find_first_intersecting_edge(mesh, intface)

Find the first interface edge crossing a mesh edge and getting inside the element

Parameters:

• mesh (2D_mesh) – mesh

• intface (disk_intface) – intface

Returns:

edge index

Return type:

int

polymemb.mesh_manager.find_next_element(elem2node, coords, cut_elems, intface, ied, cut_iel, intsec_points)

Find what is the next element cut by the intface: it is entered by the edge, it is adjacent to the cut element and it is activated

Parameters:

• elem2node (list(list(int))) – elem2node connectivity of a partially cut mesh

• coords (list(np.array)) – list of points

• intface (disk_interface) – interface

• ied (int) – idx of edge to check

• cut_iel (int) – last cut element

• intsec_points (list(int)) – list of points intersections with cut elem

Returns:

index of entered element (which is next to be cut)

Return type:

iel_to_cut (int)

polymemb.mesh_manager.get_clip(subjectPolygon, clipPolygon)

Determine clip between two polygons using Sutherland-Hodgman algorithm. (adapted from: https://rosettacode.org/wiki/Sutherland-Hodgman_polygon_clipping#Python) For reference: https://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman_algorithm

polymemb.mesh_manager.get_final_agglomerating(agglomerations, iag, verbose=False)

Internal procedure of agglomerate: Same element can find itself to be agglomerated several times, Get final element that ends up agglomerating element with index iag

Parameters:

• agglomerations (list(list(int))) – list of agglomerations to perform

• iag (int) – index in list of agglomerations

polymemb.mesh_manager.get_intersections(coords, elem, intface, ied)

Finds intersections between a simlpicial mesh element and an interface edge

Parameters:

• coords (list(np.array)) – list of points

• elem (list(np.array)) – elem to intersect (list of points)

• intface (disk_interface) – interface

• ied (int) – edge index

Returns:

intersection coordinates (if any) cut_elem_edges (list): edges that are cut

Return type:

intersections (list)

polymemb.mesh_manager.get_intface_data(mesh, iel)

Returns [cut, no_edges, first_ie] (data are used to update cuts field of agglomeratr mesh) :type mesh: mema.mesh2D :param mesh: mesh :type mesh: mema.mesh2D :type iel: int :param iel: element index :type iel: int

polymemb.mesh_manager.get_neighbor(elem2node, elem_bnd_mask, iel, v0, v1)

Look for the element sharing with iel the vertices v0 and v1

Parameters:

• elem2node (list) – current list of element, already modified with new points

• elem_bnd_mak (list) – to check if element lies on boundary

• v0 (int) – vertex 0 (index)

• v1 (int) – vertex 1 (index)

Returns:

index of found neighbor int: local index of edge to split

Return type:

int

polymemb.mesh_manager.get_triangle_vertices(elem_coords)

Given a polygon with a triangular shape (such as those that you can get after nonconformal refinement (hanging nodes)) Extracts the vertices of the triangle (local index) Beware: it only works with actual triangles, don’t apply after breaking along the mesh

polymemb.mesh_manager.glob_idx(mesh, iel, ino)

polymemb.mesh_manager.harmonize_cryterion(iel, elem_coords)

Check if element has neighbors that are too small

Parameters:

• iel (int) – element index

• elem_coords (list) – element node coordinates

Returns:

bool

polymemb.mesh_manager.initialize_disk_interface(N, rho, k_b, k_str)

Create a brand-new disk interface

Parameters:

• N (int) – number of edges/nodes

• rho (float) – radius of the disk

polymemb.mesh_manager.initialize_interface_from_stl(intface_filename, k_b, k_str)

Create the interface reading an stl file

Parameters:

• intface_filename (string) – path to stl file

• rho (float) – radius of the disk

polymemb.mesh_manager.intface_cryterion(iel, elem_coords, intface)

Check if element is intersected by interface (new version)

Parameters:

• iel (int) – element index

• elem_coords (list) – element node coordinates

• intface (mema.disk_intface) – interface

Returns:

bool

polymemb.mesh_manager.is_cut(mesh, iel, rho)

checks if face is intersected by disk with radius rho

polymemb.mesh_manager.is_proper(value)

Checks if value is real and in interval(0,1].

Parameters: - value: numerical value to check

Returns: True or False

polymemb.mesh_manager.kron(i, j)

Kronecker delta

polymemb.mesh_manager.load_square_mesh(filename, bnd_dof_type='edge', scaling=1.0)

Loads and parses a 2D mesh in msh format

Parameters: - filename: path to the mesh file - no_points_per_elem = increase size of elem2edge columns to fit polygonal elements - bnd_dof_type (logical): type of dofs on the boundary - scaling (double): scaling

Returns: - Mesh object

polymemb.mesh_manager.mark_bad_quality_elements(mesh, reference_mesh, cryt_size, cryt_skewness)

Provides a quality mask for mesh elements :type mesh: mema.2Dmesh :param mesh: mesh :type mesh: mema.2Dmesh :type reference_mesh: mema.mesh2D :param reference_mesh: uncut mesh :type reference_mesh: mema.mesh2D :type cryt_size: float :param cryt_size: minimal accepted size relative to standard far-from interface simplex :type cryt_size: float :type cryt_skewness: float :param cryt_skewness: maximal accepted aspect ratio (skewness) :type cryt_skewness: float

Returns:

0 if good quality, 1 if bad

Return type:

list (int)

class polymemb.mesh_manager.mesh2D(elem2node, coords, bnd_dof_type, d)

Bases: object

Mesh Class. Supports broken elements. It features both elem2edge and elem2node connectivity. It incorporates information about the phase and interface edges.

Attributes:

elem2node (list(list(int))): element 2 node connectivity elem2edge (list(list(int))): element 2 edge connectivity no_edges (int): number of edges edge2elem (list(list(int))): edge 2 element connectivity coords (list(np.array)): node coordinates, 2 components side_mask (list): element mask to define phase (0 int, 1 ext) node_bnd_mask (list(int)): node mask to define if on border and what border (0 if internal, k if on boundary k) edge_bnd_mask (list(int)): edge mask to define if on border elem_bnd_mask (list(int)): elem mask to define if on border cuts (list): for each cut a list in the form: [couple, starting_ie, edge2intface]

couple: elements along the cut [in, ex] starting_ie: local idx of first edge along cut[in, ex] edge2intface: map from edge on cut to idx on intface

intface_edges (list): list of interface edges (global index) d (float): semi-side of box (necessary for marking boundary nodes)

How to loop over interface edges: loop over cuts, and run over element edges using index of first and len(edge2intface)

barycenter(iel)

Determines barycenter of polygonal element (using shoelace formula)

Parameters:

iel (int) – element index

Returns:

barycenter

Return type:

np.array(float)

calc_surface(iel)

generate_edge2elem()

Calculates edge 2 element connectivity

generate_elem2edge()

Creates the local to global map for edges self.elem2edge

Parameters: - self (self2D): self

Returns: - list : elem2edge connectivity

generate_elem2elem()

Generates elem2elem connectivity (elements sharing one edge standing on same side)

generate_intface_edges()

Creates self.intface_edges (list of global indexes of interface edges) Edges are not in ascending ordered, position is correlated to order on interface (not perfectly though)

Parameters: - self (self2D): self

Returns: - list: list to replace as self.intface_edges

get_edge_length(iel, ie)

Parameters:

• iel (int) – element index

• ie (int) – index of face

Returns:

length of the edge ie

Return type:

float

Raises: ValueError: If the result is 0.

get_edge_normal(iel, num_face)

Parameters:

• self (ddr.self2D) – self

• iel (int) – element index

• num_face (int) – index of face

Returns:

edge normal

Return type:

np.array(float)

get_mesh_size()

Provides minimal and maximal element size (taken as sqr of surface)

get_xE(iel, ie)

Parameters:

• iel (int) – element index

• ie (int) – index of edge

Returns:

barycenter of edge

Return type:

float

mark_bnd_edges()

Sets mask for edges lying on boundary (no distinction between sections of boundary)

Returns: void

mark_bnd_elems()

Assigns the boundary mask for elements; It only works for square boxes

Convention:

-1: internal

1: x=-d 2: y=-d 3: x= d 4: y= d

mark_bnd_points()

Assigns the boundary mask for points; It only works for square boxes

Convention:

-1: internal

1: x=-d 2: y=-d 3: x= d 4: y= d

polymemb.mesh_manager.move_critical_points(mesh, rho)

Displace critical points which are too close to the interface (avoid rare but nasty situations of elements not properly cut)

polymemb.mesh_manager.no_check_and_add_new_points(new_coords, new_points)

Add a bunch of new points without checking if point is already there, and provide indexes on the new entries

Parameters:

• new_coords (list(np.array)) – list of points

• new_points (list(np.arra)) – points to add

polymemb.mesh_manager.point_is_inside(polygon, q)

Determine if point is inside a polygon :type polygon: list(np.array) :param polygon: polygon as list of points :type polygon: list(np.array) :type q: np.array :param q: point to check :type q: np.array

Returns:

boolean

polymemb.mesh_manager.refine_mesh(mesh0, cryterion=<function default_cryterion>, additional_arguments=[], debug=False)

Refines a triangular mesh; triangular in shape, but not necessarily as for no of points (we want to be able to realize several levels of refinement)

Parameters:

• mesh (mema.mesh2D) – mesh to refine

• cryterion (callable) – a function that specifies the cryterion to refine an element. It has signature: cryterion(iel, elem_coords, *add_args), where elem_coords is the list of coordinates of nodes of an element (as a pseudo-triangle), and add_args is a list of additional arguments

• additional_arguments (list) – addtional arguments for cryterion

• debug (bool) – print messages

Returns:

mema.mesh2D

polymemb.mesh_manager.rotate_90_clockwise(vec)

Rotate a vector 90 degrees clockwise :type vec: np.array :param vec: vector in R^2 :type vec: np.array

Returns:

rotated vector in R^2

Return type:

(np.arrya)

polymemb.mesh_manager.share_an_edge(elem2node, iel1, iel2)

Determine whether two elemnts share an edge

Parameters:

• elem2node (list(list(int))) – element2node connectivity

• iel1 (int) – first elem index

• iel2 (int) – second elem index

Returns:

boolean

polymemb.mesh_manager.shoelace_volume(coords)

Calculates volume of intface with shoelace formula (given its coords) :param coords list: formatted as intface.coords :type coords list: np.array

polymemb.mesh_manager.suppress_elements(elem2node, cuts, list_of_elems)

Suppress elements from mesh cuts field is updated accordingly :type elem2node: list(list(int)) :param elem2node: element to node connectivity :type elem2node: list(list(int)) :type cuts: list :param cuts: cuts field (interface-related data) :type cuts: list :type list_of_elems: list(int) :param list_of_elems: list of elements to suppress :type list_of_elems: list(int)

Returns:

new element to node connectivity cuts: new cuts field

Return type:

new_elem2node

polymemb.mesh_manager.visualize_intface(fig, ax, intface, color='red', show_velocity=False, show_tension=False, show_nodes=False, mark_first_node=False, show_edge_idx=False)

Plots interface

Parameters:

• fig – figure

• ax – axis

• intface (mema.disk_interface) – interface to plot

• show_velocity (boolean) – whether to show arrows for velocity

• show_tension (boolean) – whether to show arrows for tension

polymemb.mesh_manager.visualize_intface_nodal_vector(fig, ax, intface, vector_data)

Plots interface

Parameters:

• fig – figure

• ax – axis

• intface (mema.disk_interface) – interface

• vector_data – edge-wise data to plot

polymemb.mesh_manager.visualize_intface_vector(fig, ax, intface, vector_data)

Plots interface

Parameters:

• fig – figure

• ax – axis

• intface (mema.disk_interface) – interface

• vector_data – edge-wise data to plot

polymemb.mesh_manager.visualize_mesh(mesh, ax, display_no_nodes=False, display_elem_id=False, display_side=False, display_couples=False, cmap='magma', display_node_id=False, show_barycenter=False, display_nodal_dof=False)

Display mesh and related info

Parameters:

• display_no_nodes (boolean) – whther to show number of nodes of each element

• display_side (boolean) – whether to display color of side

• display_node_id (boolean) – whether to show index of nodes

• show_barycenter (boolean) – whether to show element barycenter

• display_nodal_dof (boolean) – whether to show dof numeber of points

polymemb.mesh_manager.visualize_mesh_element_data(mesh, data, ax, binary=True, cmap='magma')

Display mesh and related info

Parameters:

• data (list(float)) – element-wise scalar data

• binary (boolean) – wether the data assumes only two values

polymemb.mesh_manager.volume_conservation_postproc(vert_pos, ref_vol)

Adjust the position of vertices to restore initial volume of the surface. Displacement with respect to position of the surface after advection minimal. Solve the minimization problem:

min_{x} (x-vert_pos)**2 with constraint: volume (vert_pos)= ref_vol

For minimization, use Sequential Quadratic Programming (SQP), implemented in library scipy.optimize

Parameters:

• vert_pos (list(np.array)) – vertices formatted as field intface.coords

• ref_vol (float) – reference volume

Returns:

optimized vertices position

Return type:

list(np.array)

Module: solvers

polymemb.solvers.assemble_A(mesh, v_rec, nu_in, nu_ex)

Assembles matrix A (constant viscosity)

Parameters:

• mesh (Mesh2D) – mesh

• v_rec (list) – velocity reconstrucntion (list with elems [A,b])

• nu_in (float) – viscosity internal

• nu_ext (float) – viscosity external

Returns:

matrix A

Return type:

np.array(no_v_dofs, no_v_dofs)

Raises:

ErrorType – Description of when this error might be raised.

polymemb.solvers.assemble_A_fast(mesh, v_rec, nu_in, nu_ex)

Assembles matrix A (constant viscosity)

Parameters:

• mesh (Mesh2D) – mesh

• v_rec (list) – velocity reconstruction (list with elems [A,b])

• nu_in (float) – viscosity internal

• nu_ext (float) – viscosity external

Returns:

matrix A

Return type:

np.array(no_v_dofs, no_v_dofs)

Raises:

ErrorType – Description of when this error might be raised.

polymemb.solvers.assemble_B(mesh)

Assembles matrix B

Parameters:

mesh (ddr.Mesh) – mesh

Returns:

matrix B

Return type:

np.array(no_v_dofs, no_p_dofs)

polymemb.solvers.assemble_B_fast(mesh)

Assembles matrix B

Parameters:

mesh (ddr.Mesh) – mesh

Returns:

matrix B

Return type:

coo_matrix

polymemb.solvers.assemble_G(mesh, sigma_in, sigma_ex)

Assemble matrix G (grad-grad)

Parameters:

• mesh (mema.mesh2D) – mesh

• sigma_in (float) – diff. coefficient internal

• sigma_ex (float) – diff. coefficient external

polymemb.solvers.assemble_JP(mesh, v_rec, verbose=False)

Faster implementation of the Jump Penalization term assembly. Necessary to recover stability at lowest order (k=0). See di2020hybrid, sec. 7.6.

Parameters:

• mesh (mesh2D) – mesh

• v_rec (list) – list of velocity reconstruction format [A,b]

• verbose (bool) – verbosity flag

Returns:

JP (jump penalization matrix)

Return type:

np.array

polymemb.solvers.assemble_JP_fast(mesh, v_rec, glob_dof_table, verbose=False)

Faster implementation of the Jump Penalization term assembly. It relies on the global dof table to avoid recomputing the local to global mapping for each element. It avoids using numpy functions to speed up the process, especially in the internal loops

Parameters:

• mesh (mesh2D) – mesh

• glob_dof_table (list) – glob2loc DOF table

• v_rec (list) – list of velocity reconstruction format [A,b]

• verbose (bool) – verbosity flag

Returns:

JP (jump penalization matrix)

Return type:

np.array

polymemb.solvers.assemble_M_gamma(mesh, sigma_in, sigma_ex)

Assembles matrix M_gamma (gradient jump at interface)

Parameters:

• mesh (mema.mesh2D) – mesh

• sigma_in (float) – diff. coefficient internal

• sigma_ex (float) – diff. coefficient external

polymemb.solvers.assemble_N_gamma(mesh, sigma_in, sigma_ex)

Assembles matrix N_gamma (jump at interface of edge values)

Parameters:

• mesh (mema.mesh2D) – mesh

• sigma_in (float) – diff. coefficient internal

• sigma_ex (float) – diff. coefficient external

polymemb.solvers.assemble_S(mesh, sigma_in, sigma_ex)

Assemble matrix S (stabilization)

Parameters:

• mesh (mema.mesh2D) – mesh

• sigma_in (float) – diff. coefficient internal

• sigma_ex (float) – diff. coefficient external

polymemb.solvers.assemble_STAB(mesh)

Assembles only stabilisation component of matrix A (debugging purpose)

Parameters:

mesh (ddr.Mesh) – mesh

Returns:

matrix STAB

Return type:

np.array(no_v_dofs, no_v_dofs)

polymemb.solvers.assemble_b_J(mesh, J_datum, sigma_in, sigma_ex)

Assembles rhs B_J (jump at interface)

Parameters:

• mesh (mema.mesh2D) – mesh

• J_datum (lambda) – jump value

• sigma_in (float) – diff. coefficient internal

• sigma_ex (float) – diff. coefficient external

polymemb.solvers.assemble_b_f(mesh, f)

Assembles vector b_f associated to <f, v>

Parameters:

• mesh (ddr.Mesh) – mesh

• f (lambda function) – forcing term

Returns:

vector b

Return type:

np.array(no_v_dofs)

polymemb.solvers.assemble_b_gamma(mesh, t_gamma)

Assembles operator associated to surface tension $tau_Gamma$

Parameters:

• mesh (mema.Mesh) – mesh

• tau_gamma (list(np.array) – edge-wise surface tension

Returns:

linear system system vector tau

Return type:

np.array(no_v_dofs)

polymemb.solvers.calc_L0_error(mesh, u, ref_sol)

Calculate L0 error :type mesh: mema.mesh2D :param mesh: mesh :type mesh: mema.mesh2D :type u: np.array :param u: numeric solution :type u: np.array :type ref_sol: lambda :param ref_sol: reference solution (wrapper to be called onto (mesh, dof)) :type ref_sol: lambda

polymemb.solvers.calc_L2_error(mesh, u, ref_sol)

Calculate L2 error :type mesh: mema.mesh2D :param mesh: mesh :type mesh: mema.mesh2D :type u: np.array :param u: numeric solution :type u: np.array :type ref_sol: lambda :param ref_sol: reference solution (wrapper to be called onto (mesh, dof)) :type ref_sol: lambda

polymemb.solvers.calc_energy_error(mesh, u, ref_sol, G, N_gamma, S)

Calculate energy error :type mesh: mema.mesh2D :param mesh: mesh :type mesh: mema.mesh2D :type u: np.array :param u: numeric solution :type u: np.array :type ref_sol: lambda :param ref_sol: reference solution as function of mesh and dof :type ref_sol: lambda :type G: np.array :param G: grad.grad matrix :type G: np.array :type N_gamma: np.array :param N_gamma: interface jump penalization matrix :type N_gamma: np.array :type S: np.array :param S: stabilization matrix :type S: np.array

polymemb.solvers.calc_r2_prod_integ(v1A, v1B, v2A, v2B)

Calculates the integral over the [0,1] of the product between 2 affine functions assigned through their values at vertices 0 and 1 (first values at 0 of f1 and f2 an then values at 1) (a bit redeundant now that there is quadrature)

Parameters:

• v1A (float) – value of p1 function 1 at node A (0)

• v1B (float) – value of p1 function 1 at node B (1)

• v2A (float) – value of p1 function 2 at node A (0)

• v2B (float) – value of p1 function 2 at node A (1)

Returns:

float

polymemb.solvers.calc_t_maxwell(mesh, intface, u, eps_in, eps_ext)

Calculates Maxwell tension from solution of DDRIN :type mesh: mesh2D :param mesh: mesh :type mesh: mesh2D :type intface: disk_intface :param intface: interface :type intface: disk_intface :type u: :param u: node_wise potential :type eps_in: :param eps_in: internal permittivity :type eps_ext: :param eps_ext: external permittivity

polymemb.solvers.compute_dof_table(mesh)

Pre-compute the map loc2glob :type mesh: mesh2D :param mesh: mesh :type mesh: mesh2D

Returns:

list of local to global dof mapping for each element

Return type:

dof_table (list)

polymemb.solvers.count_dof(mesh)

Total DOF number

polymemb.solvers.disc_grad(mesh, iel)

Computes the discrete gradient matrix of an element, such that GRAD*[nodal values] returns the gradient of the local distribution of dofs [nodal values], i.e.: j-th column is discrete gradient of basis function of node j

Parameters:

• mesh (mema.mesh2D) – mesh

• iel (int) – element index

polymemb.solvers.dof_loc2glob(mesh, iel, loc_dof)

Maps local dof index to global dof index

Parameters:

• mesh (ddr.mesh2D) – mesh

• iel (int) – element index

• loc_dof (int) – local degree of freedom

Returns:

global dof index

Return type:

int

Raises:

ErrorType – Description of when this error might be raised.

polymemb.solvers.elem_velocity_energy_norm(mesh, iel, v_h)

Parameters:

• mesh (mesh2D) – mesh

• iel (int) – elem index

• v_h (np.array(float)) – velocity (coefficients in HHO space), dim=no_v_dofs (only velocity array)

polymemb.solvers.get_disc_div(mesh, iel, dof)

Provides coefficients of r^{k+1}_T for the local basis function associated to the local degree of freedom w.r.t. the canonical basis. i.e. returns [A, b] if r^{k+1}_T = A*(x-xT)+ b (x, xT, b are in R^2) In particular, A is the gradient to use in aT, and b is the average

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• dof (int) – local degree of freedom

Returns:

A (grad r^{k+1}_T) np.array (2*1): b

Return type:

np.array (2*2)

Raises:

... –

polymemb.solvers.get_disc_div_from_localV(mesh, iel, localV)

Discrete divergence

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• localV (np.array, 2*1+no_edges) – local function (v_T, {v_E}_E)

Returns:

dicrete divergence (p0(T): constant)

Return type:

float

Raises:

... –

polymemb.solvers.get_edge_dofs(mesh, iel, ie)

Description: returns 2 velocity edge dofs given element index and local edge index

Parameters:

• mesh (ddr.mesh2D) – mesh

• iel (int) – element index

• ie (int) – local index of edge

Returns:

dof_x, dof_y

Return type:

list (int)

polymemb.solvers.get_local_interpolation(mesh, iel, ref_v)

Provides local p0 polynomials associated to a velocity dof, with the convention [v_T, {v_F}_F]. v_T is p0(T)^2, v_F is p0(F)^2 for each F

Parameters:

mesh (ddr.Mesh) – mesh

iel (int): element index ref_v (lambda): H1 velocity to interpolate

Returns:

local p0 polynomials (R^2 constants) (underline{v}_T)

Return type:

np.array (2*(1+num_edge))

Raises:

... –

polymemb.solvers.get_local_polynomials(mesh, iel, dof)

Provides local p0 polynomials associated to a velocity dof, with the convention [v_T, {v_F}_F]. v_T is p0(T)^2, v_F is p0(F)^2 for each F

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• dof (int) – local degree of freedom

Returns:

local p0 polynomials (R^2 constants)

Return type:

list (np.array(2)) len=(1+num_edge))

Raises:

... –

polymemb.solvers.get_local_projections(mesh, iel, ref_v)

Given a reference velocity, computes the projection over the local basis functions and packs them into an array (usual convention for ordering). Combining get_local_polynomials and get_local_projections should be the same as get_local_interpolation

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• dof (int) – local degree of freedom

Returns:

local projections

Return type:

np.array (no_loc_dofs)

Raises:

... –

polymemb.solvers.get_local_velocities(mesh, iel, v_h)

Provides local velocities associated to element w.r.t. a global velocity unknown v_h with the convention [v_T, {v_E}_E]. v_T is p0(T)^2, v_E is p0(E)^2 for each E

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• v_h (np_array) – local degree of freedom

Returns:

local velocities [v_T, {v_E}]

Return type:

list (np.array(2)) len=(1+num_edge))

polymemb.solvers.get_v_rec(mesh, iel, dof)

Provides coefficients of r^{k+1}_T for the local basis function associated to the local degree of freedom w.r.t. the canonical basis. i.e. returns [A, b] if r^{k+1}_T = A*(x-xT)+ b (x, xT, b are in R^2) In particular, A is the gradient to use in aT, and b is the average

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• dof (int) – local degree of freedom

Returns:

A (grad r^{k+1}_T) np.array (2*1): b

Return type:

np.array (2*2)

Raises:

... –

polymemb.solvers.get_v_rec_fast(mesh, iel, dof)

Provides coefficients of r^{k+1}_T for the local basis function associated to the local degree of freedom w.r.t. the canonical basis. i.e. returns [A, b] if r^{k+1}_T = A*(x-xT)+ b (x, xT, b are in R^2) In particular, A is the gradient to use in aT, and b is the average

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• dof (int) – local degree of freedom

Returns:

A (grad r^{k+1}_T) np.array (2*1): b

Return type:

np.array (2*2)

Raises:

... –

polymemb.solvers.get_v_rec_from_localV(mesh, iel, localV)

Parameters:

• mesh (Mesh2D) – mesh

• iel (int) – element index

• list (localV) – local function (v_T, {v_E}_E)

Returns:

A (grad r^{k+1}_T) np.array (2*1): b (1/mod_T*int_T r^{k+1}_T : element average)

Return type:

np.array (2*2)

Raises:

... –

polymemb.solvers.get_v_rec_from_localV_fast(mesh, iel, localV)

Made faster by removing numpy functions.

Parameters:

• mesh (Mesh2D) – mesh

• iel (int) – element index

• list (localV) – local function (v_T, {v_E}_E)

Returns:

[A11, A12, A21, A22, B1, B2]

Return type:

tuple

polymemb.solvers.glob_idx(mesh, iel, ino)

Local DOF to global DOF connectivity for solver DDR_interfaces Convention for dof ordering: [internal unknowns, intface unknowns side in, intface unknowns side ex]

Parameters:

• mesh (mema.mesh2D) – mesh

• iel (int) – element index

• ino (int) – local node index

polymemb.solvers.impose_bc(mesh, A, b, ref_sol)

Modifies the system to impose boundary conditions

Parameters:

• mesh (mema.mesh2D) – mesh

• A (np.array) – lhs matrix

• b (np.array) – rhs array

• ref_sol (lambda) – reference solution

polymemb.solvers.impose_bc_fast(mesh, S_csr, b, ref_sol_v, zero_mean=True)

Forces Dirichlet boundary conditions on velocity Enlarges the system with one line to account for zero mean pressure

Parameters:

• mesh (ddr.Mesh) – mesh

• S (csr_matrix)

• b (np.array) – global system vector

• ref_sol_v (lambda) – reference velocity

• zero_mean (boolean) – True if you want to force zero mean for pressure

Returns: S (np.array): global system matrix

b (np.array): global system array

polymemb.solvers.interpolate_pressure(mesh, ref_p)

Parameters:

• mesh (ddr.Mesh) – mesh

• ref_p (lambda) – reference pressure

Returns:

discrete pressure

Return type:

np.array(no_p_dofs)

Remarks: L2 projection over T is replaced by value in barycenter

polymemb.solvers.interpolate_solution(mesh, ref_p, ref_v)

Parameters:

• mesh (ddr.Mesh) – mesh

• ref_p (lambda) – reference pressure

• ref_v (lambda) – reference velocity

Returns:

discrete solution in the conventinal format [p_h, v_h]

Return type:

np.array(no_tot_dofs)

Remarks: L2 projection over T/E replaced by value in barycenter

polymemb.solvers.interpolate_velocity(mesh, ref_v)

Parameters:

• mesh (ddr.Mesh) – mesh

• ref_v (lambda) – reference velocity

Returns:

discrete velocity

Return type:

np.array(no_v_dofs)

Remarks: L2 projection over T/E replaced by value in barycenter

polymemb.solvers.loc_dof_description(mesh, iel, loc_dof)

Returns description of the local dof (p/V, x/y, T/F, face_number)

Parameters:

• mesh (ddr.mesh2D) – mesh

• iel – index of elem

• loc_dof (int) – local degree of freedom

Returns:

[“p”/”V”, “-“/x”/”y”, “T”/”F”, face_number]

Return type:

list of str

Raises:

... –

polymemb.solvers.local_contribution_aT(mesh, v_rec, iel, i, j, nu)

Calculate local contribution aT (grad:grad)

Parameters:

• mesh (Mesh2D) – mesh

• v_rec (list) – list of velocity reconstrunction format [A,b]

• iel (int) – element index

• i (int) – local dof 1

• j (int) – local dof 2

• nu (float) – element viscosity

Returns:

local contribution

Return type:

real

Raises:

... –

polymemb.solvers.local_contribution_aT_fast(mesh, v_rec, iel, i, j, nu)

Calculate local contribution aT (grad:grad)

Parameters:

• mesh (Mesh2D) – mesh

• v_rec (list) – list of velocity reconstruction format tuple

• iel (int) – element index

• i (int) – local dof 1

• j (int) – local dof 2

• nu (float) – element viscosity

Returns:

local contribution

Return type:

real

Raises:

... –

polymemb.solvers.local_contribution_bT(mesh, iel, i)

Calculate local contribution bT (coupling v,p)

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• i (int) – local dof associated to v

Returns:

local contribution

Return type:

real

polymemb.solvers.local_contribution_fT(mesh, iel, f, dof)

Calculate local contribution fT (forcing term)

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• f (lambda function) – forcing term

• dof (int) – local dof

Returns:

local contribution

Return type:

real

polymemb.solvers.local_contribution_sT(mesh, v_rec, iel, i, j)

Calculate local contribution sT (stabilization)

Parameters:

• mesh (ddr.Mesh) – mesh

• v_rec (list) – list of velocity reconstrunction format [A,b]

• iel (int) – element index

• i (int) – local dof 1

• j (int) – local dof 2

Returns:

local contribution

Return type:

real

Raises:

... –

polymemb.solvers.local_contribution_sT_fast(mesh, v_rec, iel, i, j)

Calculate local contribution sT (stabilization)

Parameters:

• mesh (ddr.Mesh) – mesh

• v_rec (list) – list of velocity reconstrunction format tuple

• iel (int) – element index

• i (int) – local dof 1

• j (int) – local dof 2

Returns:

local contribution

Return type:

real

polymemb.solvers.p1_rec_matrix(mesh, iel)

Returns the P1 reconstruction matrix such that [p1 rec nodal values] = R*[nodal_values]

Parameters:

• mesh (mesh.mesh2D) – mesh

• iel (int) – element index

Returns:

np.array

polymemb.solvers.p_v_error(mesh, p_v_lambda, ref_p, ref_v, normalized=True)

Parameters:

• mesh (mesh2D) – mesh

• p_v_lambda (np.array(float)) – HHO-stokes solution (with 1 DOF for Lagrange multiplier)

• ref_p (lambda) – reference pressure

• ref_v (lambda) – reference velocity

Returns:

[pressure L2 error, velocity energy error]

Return type:

list(float)

polymemb.solvers.pressure_L2_norm(mesh, p_h)

Parameters:

• mesh (mesh2D) – mesh

• p_h (np.array(float)) – pressure (coefficients in HHO space) dim=no_p_dofs (only pressure array)

polymemb.solvers.pretty_visualize(mesh, u, fig, ax, cmap='magma')

Visualize numeric solution

Parameters:

• mesh (mema.mesh2D) – mesh

• u (np.array) – numeric solution

• fig – figure

• ax – figure

• cmap – colormap

polymemb.solvers.pretty_visualize_gradient(mesh, u, fig, ax, cmap='magma', random_threshold=3)

polymemb.solvers.reference_solution(mesh, dof, rho, ref_sol_in, ref_sol_ex, t=0, t_dep=False)

polymemb.solvers.solve(mesh, ref_sol, sigma_in, sigma_ex, eta, J_datum)

Solves an elliptic interface problem, given mesh and data.

Parameters:

• mesh (mema.mesh2D) – mesh

• ref_sol (lambda) – reference solution (or boundary data)

• sigma_in (float) – diffusion coefficient internal

• sigma_ex (float) – diffusion coefficient external

• eta (float) – user-dependent parameter for jump conditoins “a la nitsche”

• J_datum (lambda) – interface jump

Returns:

solution u np.array : global system matrix A np.array : global system r.h.s. b np.array : block matrix G (grad*grad) np.array : block matrix M_gamma (interface jump of gradient) np.array : block matrix N_gamma (interface jump) np.array : block r.h.s. b_j (interface jump)

Return type:

np.array

polymemb.solvers.solve_stokes(mesh, ref_sol_v, vol_force, nu_in, nu_ex, intface=None, with_surface_tension=False, external_tension=None, jump_penalization=False)

Solve Stokes problem with an interface and surface tension

Parameters:

• mesh (mema.Mesh2D) – mesh

• intface (mema.disk_interface) – interface

• ref_sol_v (lambda) – reference velocity

• vol_force (lambda)

• nu (float) – viscosity (up to now uniform)

• external_tension (list(np.array)) – external tension (e.g. Maxwell)

• jump_penalization (boolean) – whether o apply jump penalization

Returns:

solution (p, v, lambda) where lambda

is Lagrange multiplier

S (np.array): global system matrix (after bnd conditions) b (np.array): global system array (after bnd conditions) A (np.array): matrix grad_s:grad_s B (np.array): matrix div*q PJ(np.array): jump penalization b_gamma (np.array): global system array (after bnd conditions)

Return type:

p_v_lambda (np.array)

polymemb.solvers.solve_stokes_fast(mesh, ref_sol_v, vol_force, nu_in, nu_ex, intface=None, with_surface_tension=False, external_tension=None, jump_penalization=False, solver_type='SPSOLVE', verbose=False)

Solve Stokes problem with an interface and surface tension

Parameters:

• mesh (mema.Mesh2D) – mesh

• intface (mema.disk_interface) – interface

• ref_sol_v (lambda) – reference velocity

• vol_force (lambda)

• nu (float) – viscosity (up to now uniform)

• external_tension (list(np.array)) – external tension (e.g. Maxwell)

• jump_penalization (boolean) – whether o apply jump penalization

Returns:

solution (p, v, lambda) where lambda

is Lagrange multiplier

S (np.array): global system matrix (after bnd conditions) b (np.array): global system array (after bnd conditions) A (np.array): matrix grad_s:grad_s B (np.array): matrix div*q PJ(np.array): jump penalization b_gamma (np.array): global system array (after bnd conditions)

Return type:

p_v_lambda (np.array)

polymemb.solvers.transfer_velocity_and_advect_interface(mesh, intface, p_v_lambda, ref_vol, dt, advect=True, vol_corr=True)

Displace interface transferring velocity from mesh edges

Parameters:

• mesh (mema.Mesh2D) – mesh

• intface (mema.disk_intface) – interface

• p_v_lambda (np.array) – stokes solution format (p,v,lambda)

• dt (float) – time step

• rf_vol (float) – volume of the interface

• advect (boolean) – whether or not to advect

Returns:

moved interface

Return type:

disk_intface

polymemb.solvers.velocity_energy_norm(mesh, v_h)

Parameters:

• mesh (mesh2D) – mesh

• v_h (np.array(float)) – velocity (coefficients in HHO space), dim=no_v_dofs (only velocity array)

polymemb.solvers.visualize_solution(mesh, v_p, fig, axes, cmaps=['magma', 'viridis'], arrows='edge', arrow_density=6)

Plots solution to Stokes Problem; on first axes velocity (magnitude), on second axes pressure

Parameters:

• mesh (ddr.Mesh) – mesh

• v_p (np.array) – pressure_velocity_solution

• fig (plt.figure) – figure

• axes (plt.axes) – axes

• cmaps (colormap) – colormaps

• arrows (string) – if “edge” edge velocity is represented with arrows, otherwise element velocity.

• arrow_density (int) – from 0 to 6 to increase arrows that are displayed

• fig – if given

• ax – if given

Returns: void

HHO Stokes

polymemb.solvers.hho.assemble_A(mesh, v_rec, nu_in, nu_ex)

Assembles matrix A (constant viscosity)

Parameters:

• mesh (Mesh2D) – mesh

• v_rec (list) – velocity reconstrucntion (list with elems [A,b])

• nu_in (float) – viscosity internal

• nu_ext (float) – viscosity external

Returns:

matrix A

Return type:

np.array(no_v_dofs, no_v_dofs)

Raises:

ErrorType – Description of when this error might be raised.

polymemb.solvers.hho.assemble_A_fast(mesh, v_rec, nu_in, nu_ex)

Assembles matrix A (constant viscosity)

Parameters:

• mesh (Mesh2D) – mesh

• v_rec (list) – velocity reconstruction (list with elems [A,b])

• nu_in (float) – viscosity internal

• nu_ext (float) – viscosity external

Returns:

matrix A

Return type:

np.array(no_v_dofs, no_v_dofs)

Raises:

ErrorType – Description of when this error might be raised.

polymemb.solvers.hho.assemble_B(mesh)

Assembles matrix B

Parameters:

mesh (ddr.Mesh) – mesh

Returns:

matrix B

Return type:

np.array(no_v_dofs, no_p_dofs)

polymemb.solvers.hho.assemble_B_fast(mesh)

Assembles matrix B

Parameters:

mesh (ddr.Mesh) – mesh

Returns:

matrix B

Return type:

coo_matrix

polymemb.solvers.hho.assemble_JP(mesh, v_rec, verbose=False)

Faster implementation of the Jump Penalization term assembly. Necessary to recover stability at lowest order (k=0). See di2020hybrid, sec. 7.6.

Parameters:

• mesh (mesh2D) – mesh

• v_rec (list) – list of velocity reconstruction format [A,b]

• verbose (bool) – verbosity flag

Returns:

JP (jump penalization matrix)

Return type:

np.array

polymemb.solvers.hho.assemble_JP_fast(mesh, v_rec, glob_dof_table, verbose=False)

Faster implementation of the Jump Penalization term assembly. It relies on the global dof table to avoid recomputing the local to global mapping for each element. It avoids using numpy functions to speed up the process, especially in the internal loops

Parameters:

• mesh (mesh2D) – mesh

• glob_dof_table (list) – glob2loc DOF table

• v_rec (list) – list of velocity reconstruction format [A,b]

• verbose (bool) – verbosity flag

Returns:

JP (jump penalization matrix)

Return type:

np.array

polymemb.solvers.hho.assemble_STAB(mesh)

Assembles only stabilisation component of matrix A (debugging purpose)

Parameters:

mesh (ddr.Mesh) – mesh

Returns:

matrix STAB

Return type:

np.array(no_v_dofs, no_v_dofs)

polymemb.solvers.hho.assemble_b_f(mesh, f)

Assembles vector b_f associated to <f, v>

Parameters:

• mesh (ddr.Mesh) – mesh

• f (lambda function) – forcing term

Returns:

vector b

Return type:

np.array(no_v_dofs)

polymemb.solvers.hho.assemble_b_gamma(mesh, t_gamma)

Assembles operator associated to surface tension $tau_Gamma$

Parameters:

• mesh (mema.Mesh) – mesh

• tau_gamma (list(np.array) – edge-wise surface tension

Returns:

linear system system vector tau

Return type:

np.array(no_v_dofs)

polymemb.solvers.hho.compute_dof_table(mesh)

Pre-compute the map loc2glob :type mesh: mesh2D :param mesh: mesh :type mesh: mesh2D

Returns:

list of local to global dof mapping for each element

Return type:

dof_table (list)

polymemb.solvers.hho.count_dof(mesh)

Calculate total number of degrees of freedom of the solver (size of the linear system)

Parameters:

mesh (ddr.mesh2D) – mesh

Returns:

number of dofs for pressure, for velocity and linear system size

Return type:

list (int)

Raises:

ErrorType – Description of when this error might be raised.

polymemb.solvers.hho.dof_loc2glob(mesh, iel, loc_dof)

Maps local dof index to global dof index

Parameters:

• mesh (ddr.mesh2D) – mesh

• iel (int) – element index

• loc_dof (int) – local degree of freedom

Returns:

global dof index

Return type:

int

Raises:

ErrorType – Description of when this error might be raised.

polymemb.solvers.hho.elem_velocity_energy_norm(mesh, iel, v_h)

Parameters:

• mesh (mesh2D) – mesh

• iel (int) – elem index

• v_h (np.array(float)) – velocity (coefficients in HHO space), dim=no_v_dofs (only velocity array)

polymemb.solvers.hho.get_disc_div(mesh, iel, dof)

Provides coefficients of r^{k+1}_T for the local basis function associated to the local degree of freedom w.r.t. the canonical basis. i.e. returns [A, b] if r^{k+1}_T = A*(x-xT)+ b (x, xT, b are in R^2) In particular, A is the gradient to use in aT, and b is the average

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• dof (int) – local degree of freedom

Returns:

A (grad r^{k+1}_T) np.array (2*1): b

Return type:

np.array (2*2)

Raises:

... –

polymemb.solvers.hho.get_disc_div_from_localV(mesh, iel, localV)

Discrete divergence

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• localV (np.array, 2*1+no_edges) – local function (v_T, {v_E}_E)

Returns:

dicrete divergence (p0(T): constant)

Return type:

float

Raises:

... –

polymemb.solvers.hho.get_edge_dofs(mesh, iel, ie)

Description: returns 2 velocity edge dofs given element index and local edge index

Parameters:

• mesh (ddr.mesh2D) – mesh

• iel (int) – element index

• ie (int) – local index of edge

Returns:

dof_x, dof_y

Return type:

list (int)

polymemb.solvers.hho.get_local_interpolation(mesh, iel, ref_v)

Provides local p0 polynomials associated to a velocity dof, with the convention [v_T, {v_F}_F]. v_T is p0(T)^2, v_F is p0(F)^2 for each F

Parameters:

mesh (ddr.Mesh) – mesh

iel (int): element index ref_v (lambda): H1 velocity to interpolate

Returns:

local p0 polynomials (R^2 constants) (underline{v}_T)

Return type:

np.array (2*(1+num_edge))

Raises:

... –

polymemb.solvers.hho.get_local_polynomials(mesh, iel, dof)

Provides local p0 polynomials associated to a velocity dof, with the convention [v_T, {v_F}_F]. v_T is p0(T)^2, v_F is p0(F)^2 for each F

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• dof (int) – local degree of freedom

Returns:

local p0 polynomials (R^2 constants)

Return type:

list (np.array(2)) len=(1+num_edge))

Raises:

... –

polymemb.solvers.hho.get_local_projections(mesh, iel, ref_v)

Given a reference velocity, computes the projection over the local basis functions and packs them into an array (usual convention for ordering). Combining get_local_polynomials and get_local_projections should be the same as get_local_interpolation

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• dof (int) – local degree of freedom

Returns:

local projections

Return type:

np.array (no_loc_dofs)

Raises:

... –

polymemb.solvers.hho.get_local_velocities(mesh, iel, v_h)

Provides local velocities associated to element w.r.t. a global velocity unknown v_h with the convention [v_T, {v_E}_E]. v_T is p0(T)^2, v_E is p0(E)^2 for each E

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• v_h (np_array) – local degree of freedom

Returns:

local velocities [v_T, {v_E}]

Return type:

list (np.array(2)) len=(1+num_edge))

polymemb.solvers.hho.get_v_rec(mesh, iel, dof)

Provides coefficients of r^{k+1}_T for the local basis function associated to the local degree of freedom w.r.t. the canonical basis. i.e. returns [A, b] if r^{k+1}_T = A*(x-xT)+ b (x, xT, b are in R^2) In particular, A is the gradient to use in aT, and b is the average

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• dof (int) – local degree of freedom

Returns:

A (grad r^{k+1}_T) np.array (2*1): b

Return type:

np.array (2*2)

Raises:

... –

polymemb.solvers.hho.get_v_rec_fast(mesh, iel, dof)

Provides coefficients of r^{k+1}_T for the local basis function associated to the local degree of freedom w.r.t. the canonical basis. i.e. returns [A, b] if r^{k+1}_T = A*(x-xT)+ b (x, xT, b are in R^2) In particular, A is the gradient to use in aT, and b is the average

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• dof (int) – local degree of freedom

Returns:

A (grad r^{k+1}_T) np.array (2*1): b

Return type:

np.array (2*2)

Raises:

... –

polymemb.solvers.hho.get_v_rec_from_localV(mesh, iel, localV)

Parameters:

• mesh (Mesh2D) – mesh

• iel (int) – element index

• list (localV) – local function (v_T, {v_E}_E)

Returns:

A (grad r^{k+1}_T) np.array (2*1): b (1/mod_T*int_T r^{k+1}_T : element average)

Return type:

np.array (2*2)

Raises:

... –

polymemb.solvers.hho.get_v_rec_from_localV_fast(mesh, iel, localV)

Made faster by removing numpy functions.

Parameters:

• mesh (Mesh2D) – mesh

• iel (int) – element index

• list (localV) – local function (v_T, {v_E}_E)

Returns:

[A11, A12, A21, A22, B1, B2]

Return type:

tuple

polymemb.solvers.hho.impose_bc(mesh, S, b, ref_sol_v, zero_mean=True)

Forces Dirichlet boundary conditions on velocity Enlarges the system with one line to account for zero mean pressure

Parameters:

• mesh (ddr.Mesh) – mesh

• S (np.array) – global system matrix

• b (np.array) – global system vector

• ref_sol_v (lambda) – reference velocity

• zero_mean (boolean) – True if you want to force zero mean for pressure

Returns: S (np.array): global system matrix

b (np.array): global system array

polymemb.solvers.hho.impose_bc_fast(mesh, S_csr, b, ref_sol_v, zero_mean=True)

Forces Dirichlet boundary conditions on velocity Enlarges the system with one line to account for zero mean pressure

Parameters:

• mesh (ddr.Mesh) – mesh

• S (csr_matrix)

• b (np.array) – global system vector

• ref_sol_v (lambda) – reference velocity

• zero_mean (boolean) – True if you want to force zero mean for pressure

Returns: S (np.array): global system matrix

b (np.array): global system array

polymemb.solvers.hho.interpolate_pressure(mesh, ref_p)

Parameters:

• mesh (ddr.Mesh) – mesh

• ref_p (lambda) – reference pressure

Returns:

discrete pressure

Return type:

np.array(no_p_dofs)

Remarks: L2 projection over T is replaced by value in barycenter

polymemb.solvers.hho.interpolate_solution(mesh, ref_p, ref_v)

Parameters:

• mesh (ddr.Mesh) – mesh

• ref_p (lambda) – reference pressure

• ref_v (lambda) – reference velocity

Returns:

discrete solution in the conventinal format [p_h, v_h]

Return type:

np.array(no_tot_dofs)

Remarks: L2 projection over T/E replaced by value in barycenter

polymemb.solvers.hho.interpolate_velocity(mesh, ref_v)

Parameters:

• mesh (ddr.Mesh) – mesh

• ref_v (lambda) – reference velocity

Returns:

discrete velocity

Return type:

np.array(no_v_dofs)

Remarks: L2 projection over T/E replaced by value in barycenter

polymemb.solvers.hho.loc_dof_description(mesh, iel, loc_dof)

Returns description of the local dof (p/V, x/y, T/F, face_number)

Parameters:

• mesh (ddr.mesh2D) – mesh

• iel – index of elem

• loc_dof (int) – local degree of freedom

Returns:

[“p”/”V”, “-“/x”/”y”, “T”/”F”, face_number]

Return type:

list of str

Raises:

... –

polymemb.solvers.hho.local_contribution_aT(mesh, v_rec, iel, i, j, nu)

Calculate local contribution aT (grad:grad)

Parameters:

• mesh (Mesh2D) – mesh

• v_rec (list) – list of velocity reconstrunction format [A,b]

• iel (int) – element index

• i (int) – local dof 1

• j (int) – local dof 2

• nu (float) – element viscosity

Returns:

local contribution

Return type:

real

Raises:

... –

polymemb.solvers.hho.local_contribution_aT_fast(mesh, v_rec, iel, i, j, nu)

Calculate local contribution aT (grad:grad)

Parameters:

• mesh (Mesh2D) – mesh

• v_rec (list) – list of velocity reconstruction format tuple

• iel (int) – element index

• i (int) – local dof 1

• j (int) – local dof 2

• nu (float) – element viscosity

Returns:

local contribution

Return type:

real

Raises:

... –

polymemb.solvers.hho.local_contribution_bT(mesh, iel, i)

Calculate local contribution bT (coupling v,p)

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• i (int) – local dof associated to v

Returns:

local contribution

Return type:

real

polymemb.solvers.hho.local_contribution_fT(mesh, iel, f, dof)

Calculate local contribution fT (forcing term)

Parameters:

• mesh (ddr.Mesh) – mesh

• iel (int) – element index

• f (lambda function) – forcing term

• dof (int) – local dof

Returns:

local contribution

Return type:

real

polymemb.solvers.hho.local_contribution_sT(mesh, v_rec, iel, i, j)

Calculate local contribution sT (stabilization)

Parameters:

• mesh (ddr.Mesh) – mesh

• v_rec (list) – list of velocity reconstrunction format [A,b]

• iel (int) – element index

• i (int) – local dof 1

• j (int) – local dof 2

Returns:

local contribution

Return type:

real

Raises:

... –

polymemb.solvers.hho.local_contribution_sT_fast(mesh, v_rec, iel, i, j)

Calculate local contribution sT (stabilization)

Parameters:

• mesh (ddr.Mesh) – mesh

• v_rec (list) – list of velocity reconstrunction format tuple

• iel (int) – element index

• i (int) – local dof 1

• j (int) – local dof 2

Returns:

local contribution

Return type:

real

polymemb.solvers.hho.p_v_error(mesh, p_v_lambda, ref_p, ref_v, normalized=True)

Parameters:

• mesh (mesh2D) – mesh

• p_v_lambda (np.array(float)) – HHO-stokes solution (with 1 DOF for Lagrange multiplier)

• ref_p (lambda) – reference pressure

• ref_v (lambda) – reference velocity

Returns:

[pressure L2 error, velocity energy error]

Return type:

list(float)

polymemb.solvers.hho.pressure_L2_norm(mesh, p_h)

Parameters:

• mesh (mesh2D) – mesh

• p_h (np.array(float)) – pressure (coefficients in HHO space) dim=no_p_dofs (only pressure array)

polymemb.solvers.hho.solve_stokes(mesh, ref_sol_v, vol_force, nu_in, nu_ex, intface=None, with_surface_tension=False, external_tension=None, jump_penalization=False)

Solve Stokes problem with an interface and surface tension

Parameters:

• mesh (mema.Mesh2D) – mesh

• intface (mema.disk_interface) – interface

• ref_sol_v (lambda) – reference velocity

• vol_force (lambda)

• nu (float) – viscosity (up to now uniform)

• external_tension (list(np.array)) – external tension (e.g. Maxwell)

• jump_penalization (boolean) – whether o apply jump penalization

Returns:

solution (p, v, lambda) where lambda

is Lagrange multiplier

S (np.array): global system matrix (after bnd conditions) b (np.array): global system array (after bnd conditions) A (np.array): matrix grad_s:grad_s B (np.array): matrix div*q PJ(np.array): jump penalization b_gamma (np.array): global system array (after bnd conditions)

Return type:

p_v_lambda (np.array)

polymemb.solvers.hho.solve_stokes_fast(mesh, ref_sol_v, vol_force, nu_in, nu_ex, intface=None, with_surface_tension=False, external_tension=None, jump_penalization=False, solver_type='SPSOLVE', verbose=False)

Solve Stokes problem with an interface and surface tension

Parameters:

• mesh (mema.Mesh2D) – mesh

• intface (mema.disk_interface) – interface

• ref_sol_v (lambda) – reference velocity

• vol_force (lambda)

• nu (float) – viscosity (up to now uniform)

• external_tension (list(np.array)) – external tension (e.g. Maxwell)

• jump_penalization (boolean) – whether o apply jump penalization

Returns:

solution (p, v, lambda) where lambda

is Lagrange multiplier

S (np.array): global system matrix (after bnd conditions) b (np.array): global system array (after bnd conditions) A (np.array): matrix grad_s:grad_s B (np.array): matrix div*q PJ(np.array): jump penalization b_gamma (np.array): global system array (after bnd conditions)

Return type:

p_v_lambda (np.array)

polymemb.solvers.hho.transfer_velocity_and_advect_interface(mesh, intface, p_v_lambda, ref_vol, dt, advect=True, vol_corr=True)

Displace interface transferring velocity from mesh edges

Parameters:

• mesh (mema.Mesh2D) – mesh

• intface (mema.disk_intface) – interface

• p_v_lambda (np.array) – stokes solution format (p,v,lambda)

• dt (float) – time step

• rf_vol (float) – volume of the interface

• advect (boolean) – whether or not to advect

Returns:

moved interface

Return type:

disk_intface

polymemb.solvers.hho.velocity_energy_norm(mesh, v_h)

Parameters:

• mesh (mesh2D) – mesh

• v_h (np.array(float)) – velocity (coefficients in HHO space), dim=no_v_dofs (only velocity array)

polymemb.solvers.hho.visualize_solution(mesh, v_p, fig, axes, cmaps=['magma', 'viridis'], arrows='edge', arrow_density=6)

Plots solution to Stokes Problem; on first axes velocity (magnitude), on second axes pressure

Parameters:

• mesh (ddr.Mesh) – mesh

• v_p (np.array) – pressure_velocity_solution

• fig (plt.figure) – figure

• axes (plt.axes) – axes

• cmaps (colormap) – colormaps

• arrows (string) – if “edge” edge velocity is represented with arrows, otherwise element velocity.

• arrow_density (int) – from 0 to 6 to increase arrows that are displayed

• fig – if given

• ax – if given

Returns: void

DDR Elliptic Interface

polymemb.solvers.ddr.assemble_G(mesh, sigma_in, sigma_ex)

Assemble matrix G (grad-grad)

Parameters:

• mesh (mema.mesh2D) – mesh

• sigma_in (float) – diff. coefficient internal

• sigma_ex (float) – diff. coefficient external

polymemb.solvers.ddr.assemble_M_gamma(mesh, sigma_in, sigma_ex)

Assembles matrix M_gamma (gradient jump at interface)

Parameters:

• mesh (mema.mesh2D) – mesh

• sigma_in (float) – diff. coefficient internal

• sigma_ex (float) – diff. coefficient external

polymemb.solvers.ddr.assemble_N_gamma(mesh, sigma_in, sigma_ex)

Assembles matrix N_gamma (jump at interface of edge values)

Parameters:

• mesh (mema.mesh2D) – mesh

• sigma_in (float) – diff. coefficient internal

• sigma_ex (float) – diff. coefficient external

polymemb.solvers.ddr.assemble_S(mesh, sigma_in, sigma_ex)

Assemble matrix S (stabilization)

Parameters:

• mesh (mema.mesh2D) – mesh

• sigma_in (float) – diff. coefficient internal

• sigma_ex (float) – diff. coefficient external

polymemb.solvers.ddr.assemble_b_J(mesh, J_datum, sigma_in, sigma_ex)

Assembles rhs B_J (jump at interface)

Parameters:

• mesh (mema.mesh2D) – mesh

• J_datum (lambda) – jump value

• sigma_in (float) – diff. coefficient internal

• sigma_ex (float) – diff. coefficient external

polymemb.solvers.ddr.calc_L0_error(mesh, u, ref_sol)

Calculate L0 error :type mesh: mema.mesh2D :param mesh: mesh :type mesh: mema.mesh2D :type u: np.array :param u: numeric solution :type u: np.array :type ref_sol: lambda :param ref_sol: reference solution (wrapper to be called onto (mesh, dof)) :type ref_sol: lambda

polymemb.solvers.ddr.calc_L2_error(mesh, u, ref_sol)

Calculate L2 error :type mesh: mema.mesh2D :param mesh: mesh :type mesh: mema.mesh2D :type u: np.array :param u: numeric solution :type u: np.array :type ref_sol: lambda :param ref_sol: reference solution (wrapper to be called onto (mesh, dof)) :type ref_sol: lambda

polymemb.solvers.ddr.calc_energy_error(mesh, u, ref_sol, G, N_gamma, S)

Calculate energy error :type mesh: mema.mesh2D :param mesh: mesh :type mesh: mema.mesh2D :type u: np.array :param u: numeric solution :type u: np.array :type ref_sol: lambda :param ref_sol: reference solution as function of mesh and dof :type ref_sol: lambda :type G: np.array :param G: grad.grad matrix :type G: np.array :type N_gamma: np.array :param N_gamma: interface jump penalization matrix :type N_gamma: np.array :type S: np.array :param S: stabilization matrix :type S: np.array

polymemb.solvers.ddr.calc_r2_prod_integ(v1A, v1B, v2A, v2B)

Calculates the integral over the [0,1] of the product between 2 affine functions assigned through their values at vertices 0 and 1 (first values at 0 of f1 and f2 an then values at 1) (a bit redeundant now that there is quadrature)

Parameters:

• v1A (float) – value of p1 function 1 at node A (0)

• v1B (float) – value of p1 function 1 at node B (1)

• v2A (float) – value of p1 function 2 at node A (0)

• v2B (float) – value of p1 function 2 at node A (1)

Returns:

float

polymemb.solvers.ddr.calc_t_maxwell(mesh, intface, u, eps_in, eps_ext)

Calculates Maxwell tension from solution of DDRIN :type mesh: mesh2D :param mesh: mesh :type mesh: mesh2D :type intface: disk_intface :param intface: interface :type intface: disk_intface :type u: :param u: node_wise potential :type eps_in: :param eps_in: internal permittivity :type eps_ext: :param eps_ext: external permittivity

polymemb.solvers.ddr.count_dof(mesh)

Total DOF number

polymemb.solvers.ddr.disc_grad(mesh, iel)

Computes the discrete gradient matrix of an element, such that GRAD*[nodal values] returns the gradient of the local distribution of dofs [nodal values], i.e.: j-th column is discrete gradient of basis function of node j

Parameters:

• mesh (mema.mesh2D) – mesh

• iel (int) – element index

polymemb.solvers.ddr.glob_idx(mesh, iel, ino)

Local DOF to global DOF connectivity for solver DDR_interfaces Convention for dof ordering: [internal unknowns, intface unknowns side in, intface unknowns side ex]

Parameters:

• mesh (mema.mesh2D) – mesh

• iel (int) – element index

• ino (int) – local node index

polymemb.solvers.ddr.impose_bc(mesh, A, b, ref_sol)

Modifies the system to impose boundary conditions

Parameters:

• mesh (mema.mesh2D) – mesh

• A (np.array) – lhs matrix

• b (np.array) – rhs array

• ref_sol (lambda) – reference solution

polymemb.solvers.ddr.p1_rec_matrix(mesh, iel)

Returns the P1 reconstruction matrix such that [p1 rec nodal values] = R*[nodal_values]

Parameters:

• mesh (mesh.mesh2D) – mesh

• iel (int) – element index

Returns:

np.array

polymemb.solvers.ddr.pretty_visualize(mesh, u, fig, ax, cmap='magma')

Visualize numeric solution

Parameters:

• mesh (mema.mesh2D) – mesh

• u (np.array) – numeric solution

• fig – figure

• ax – figure

• cmap – colormap

polymemb.solvers.ddr.pretty_visualize_gradient(mesh, u, fig, ax, cmap='magma', random_threshold=3)

polymemb.solvers.ddr.reference_solution(mesh, dof, rho, ref_sol_in, ref_sol_ex, t=0, t_dep=False)

polymemb.solvers.ddr.solve(mesh, ref_sol, sigma_in, sigma_ex, eta, J_datum)

Solves an elliptic interface problem, given mesh and data.

Parameters:

• mesh (mema.mesh2D) – mesh

• ref_sol (lambda) – reference solution (or boundary data)

• sigma_in (float) – diffusion coefficient internal

• sigma_ex (float) – diffusion coefficient external

• eta (float) – user-dependent parameter for jump conditoins “a la nitsche”

• J_datum (lambda) – interface jump

Returns:

solution u np.array : global system matrix A np.array : global system r.h.s. b np.array : block matrix G (grad*grad) np.array : block matrix M_gamma (interface jump of gradient) np.array : block matrix N_gamma (interface jump) np.array : block r.h.s. b_j (interface jump)

Return type:

np.array

Bibliography

[1]

D. A. Di Pietro and J. Droniou. The Hybrid High-Order method for polytopal meshes. Number 19 in Modeling, Simulation and Application. Springer International Publishing, 2020. doi:10.1007/978-3-030-37203-3.

[2]

Daniele A. Di Pietro and Jérôme Droniou. An arbitrary-order discrete de Rham complex on polyhedral meshes: Exactness, Poincaré inequalities, and consistency. Found. Comput. Math., 23:85–164, 2023. doi:10.1007/s10208-021-09542-8.

Indices and tables

• Index

• Module Index

• Search Page