Different results with MPI

General
1

Hello,
I have problems in understanding how mpi works. In the code below, I am trying to calculate the error of FEA results with the exact solution. Also I want to plot the displacements of a certain node through time. But when I run in parallel I get different results:

import dolfinx as fe
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem, assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import os
from mpi4py import MPI
import ufl
from petsc4py import PETSc
import pyvista
import pandas as pd
import meshio
from scipy.spatial import cKDTree
from scipy.special import jv, yv, iv, kv
from itertools import product
import h5py
import matplotlib.pyplot as plt


def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data.astype(np.int32)]})
    return out_mesh


def main(rho, E, nu, maximum_pressure, minimum_pressure, mesh_file, num_steps=100, element_degree=1):

    maximum_pressure = fe.default_scalar_type(maximum_pressure)
    minimum_pressure = fe.default_scalar_type(minimum_pressure)
    rho = fe.default_scalar_type(rho)
    E = fe.default_scalar_type(E)
    nu = fe.default_scalar_type(nu)

    proc = MPI.COMM_WORLD.rank

    ############################################Importing the mesh############################################
    ##########################################################################################################
    #Wrinting the mesh as XDMF file to read from that later
    if proc == 0:
        # Read in mesh
        msh = meshio.read(mesh_file)

        # Create and save one file for the mesh, and one file for the facets
        tetra_mesh = create_mesh(msh, "tetra", prune_z=False)
        triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
        facets = create_mesh(msh, "triangle", prune_z=False)
        meshio.write(os.path.join("output", "mt.xdmf"), triangle_mesh)
        meshio.write(os.path.join("output", "mesh.xdmf"), tetra_mesh)
    
    MPI.COMM_WORLD.barrier()

    #Reading the xdmf files
    with fe.io.XDMFFile(MPI.COMM_WORLD, os.path.join("output", "mesh.xdmf"), "r") as xdmf:
        mesh = xdmf.read_mesh(name="Grid")
        cell_markers = xdmf.read_meshtags(mesh, name="Grid")

    gdim = mesh.topology.dim
    mesh.topology.create_connectivity(gdim, gdim - 1)
    with fe.io.XDMFFile(MPI.COMM_WORLD, os.path.join("output", "mt.xdmf"), "r") as xdmf:
        facet_markers = xdmf.read_meshtags(mesh, name="Grid")

    points = mesh.geometry.x

    radius_to_check = 2
    for i, p in enumerate(points):
        if np.linalg.norm(p) == radius_to_check:
            index_of_radius_to_check = i
            break

    ######################################################################################################
    #######################creating function space and boundary condition#################################
    element = ufl.VectorElement("Lagrange", mesh.ufl_cell(), element_degree)
    V = fe.fem.functionspace(mesh, element)

    gamma_p = 0.01
    t = np.linspace(0, 2*np.pi / gamma_p, num_steps)
    dt = fe.fem.Constant(mesh, t[1] - t[0])

    pressure = maximum_pressure / 2 * (1 - np.cos(gamma_p * t))
    p = fe.fem.Constant(mesh, pressure[0])   

    ######################################################################################################
    ######################################################################################################

    ######################################################################################################
    #########################################Formulation##################################################
    u = ufl.TrialFunction(V)           
    v  = ufl.TestFunction(V)          
    old_u  = fe.fem.Function(V)
    old_velocity  = fe.fem.Function(V)
    old_acceleration  = fe.fem.Function(V)

    d = len(u)
    I = ufl.variable(ufl.Identity(d))             # Identity tensor
    F = ufl.variable(I + ufl.grad(u))             # Deformation gradient
    C = ufl.variable(F.T*F)                   # Right Cauchy-Green tensor
    J = ufl.variable(ufl.det(F))

    metadata = {"quadrature_degree": 2}   
    ds = ufl.Measure('ds', domain=mesh, subdomain_data=facet_markers, metadata=metadata)
    dx = ufl.Measure("dx", domain=mesh, metadata=metadata)

    # Generalized-alpha method parameters
    alpha_m = fe.fem.Constant(mesh, 0.01) 
    alpha_f = fe.fem.Constant(mesh, 0.01) 
    gamma   = fe.fem.Constant(mesh, fe.default_scalar_type(0.5+alpha_f-alpha_m))
    beta    = fe.fem.Constant(mesh, fe.default_scalar_type((gamma+0.5)**2/4.))

    # Update formula for acceleration
    # a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0)
    def update_a(u, u_old, v_old, a_old, ufl=True):
        if ufl:
            dt_ = dt
            beta_ = beta
        else:
            dt_ = float(dt)
            beta_ = float(beta)
        return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old

    # Update formula for velocity
    # v = dt * ((1-gamma)*a0 + gamma*a) + v0
    def update_v(a, u_old, v_old, a_old, ufl=True):
        if ufl:
            dt_ = dt
            gamma_ = gamma
        else:
            dt_ = float(dt)
            gamma_ = float(gamma)
        return v_old + dt_*((1-gamma_)*a_old + gamma_*a)


    def update_fields(u_new, u_old, v_old, a_old):
        """Update fields at the end of each time step."""
        # Get vectors (references)
        u_vec, u0_vec  = u_new.x.array[:], u_old.x.array[:]
        v0_vec, a0_vec = v_old.x.array[:], a_old.x.array[:]

        # use update functions using vector arguments
        a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False)
        v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False)

        # Update (u_old <- u)
        v_old.x.array[:], a_old.x.array[:] = v_vec, a_vec
        u_old.x.array[:] = u_new.x.array[:]


    def avg(x_old, x_new, alpha):
        return alpha*x_old + (1-alpha)*x_new

    normal = -ufl.FacetNormal(mesh)

    mu = fe.fem.Constant(mesh, E / (2 * (1 + nu)))
    lmbda = fe.fem.Constant(mesh, E * nu / ((1 + nu) * (1 - 2 * nu)))
    
    def epsilon(u):
        return ufl.sym(ufl.grad(u))

    def S(u):
        return 2.0 * mu * epsilon(u) + lmbda * ufl.nabla_div(u) * I 

    acceleration = update_a(u, old_u, old_velocity, old_acceleration, ufl=False)  
    velocity = update_v(acceleration, old_u, old_velocity, old_acceleration, ufl=True)

    formulation = rho * ufl.inner(avg(old_acceleration, acceleration, alpha_m), v) * dx \
      + ufl.inner(epsilon(v), S(avg(old_u, u, alpha_f))) * dx \
      - ufl.dot(v, p * normal) * ds(112)   

    bilinear_form = fe.fem.form(ufl.lhs(formulation))
    linear_form = fe.fem.form(ufl.rhs(formulation))

    ######################################################################################################
    ######################################################################################################

    ######################################################################################################
    ###############################################Solving################################################
    A = assemble_matrix(bilinear_form, bcs=[])
    A.assemble()
    b = create_vector(linear_form)

    solver = PETSc.KSP().create(mesh.comm)
    solver.setInitialGuessNonzero(True)
    solver.setOperators(A)
    solver.getPC().setType(PETSc.PC.Type.SOR)

    displacement = fe.fem.Function(V)

    class exact_solution():
        def __init__(self, rho, inner_radius, k, lmbda, miu, p_0, gamma, num_steps):  # p = -p_0 * (1-cos(gammma*tau))
            self.a = inner_radius
            self.k = k  # outer_radius = k * a
            self.phi = -1/2
            self.n = 3/2
            self.C_11 = lmbda + 2 * miu #E * (1 - miu) / ((1 + miu) * (1 - 2*miu))
            self.C_12 = lmbda #E * miu / ((1 + miu) * (1 - 2*miu))
            self.m = self.C_12 / self.C_11
            self.S_1 = self.n - 2 * self.m - self.phi
            self.alpha = [1.01034, 3.33252, 9.48585, 6.37556, 15.74444, 12.61203, 28.29456, 18.87993, 25.1555, 22.01717, 31.43413, 34.57406, 37.71428, 40.8547, 43.99529, 47.13602, 50.27685, 62.84095, 53.41778, 59.69984, 56.55877, 65.98211, 78.54709, 69.12331, 109.96094, 72.26454, 75.4058, 84.82974, 87.97109, 81.68841, 144.51722, 97.39524, 94.25384, 91.11246, 141.37571, 100.53665, 103.67807, 106.8195, 113.10239, 116.24384, 128.80974, 119.38531]
            self.alpha.sort()
            self.p_0 = p_0
            self.gamma = gamma
            self.tau = 0 # 6 * 2 * np.pi
            self.c = np.sqrt(self.C_11 / rho)
            self.dtau =  2*np.pi / self.gamma /num_steps # dt /(2*np.pi)  # * self.c


        def return_tau(self):
            return self.tau


        def step_one_period(self):
            self.tau += self.gamma * 2*np.pi # self.c 


        def step_in_time(self):
            self.tau += self.dtau
            # print(f"tau: {self.tau}")


        def r(self, x):   # distance to the center 
            return np.sqrt(np.power(x[0], 2) + np.power(x[1], 2) + np.power(x[2], 2))


        def theta_cood(self, x):
            return np.arccos(x[2]/self.r(x))


        def phi_cood(self, x):
            return np.sign(x[1]) * np.arccos(x[0]/np.sqrt(np.power(x[0], 2) + np.power(x[1], 2)))


        def F(self, p, x):
            k = self.k
            n = self.n
            S_1 = self.S_1
            a = self.a
            return (kv(n, p*k + 0j) * S_1 + p*k * kv(n-1, p*k + 0j)) * iv(n, p*self.r(x)/a + 0j)\
                   - (iv(n, p*k + 0j) * S_1 - p*k * iv(n-1, p*k + 0j)) * kv(n, p*self.r(x)/a + 0j)


        def G(self, p):
            k = self.k
            n = self.n
            S_1 = self.S_1
            a = self.a
            return (kv(n, p*k + 0j) * S_1 + p*k * kv(n-1, p*k + 0j)) * (iv(n, p + 0j) * S_1 - p * iv(n-1, p + 0j))\
                   - (kv(n, p + 0j) * S_1 + p * kv(n-1, p + 0j)) * (iv(n, p*k + 0j) * S_1 - p*k * iv(n-1, p*k + 0j))


        def __call__(self, x, radial=False):

            p_0 = self.p_0
            a = self.a
            phi = self.phi
            C_11 = self.C_11
            k = self.k
            n = self.n
            S_1 = self.S_1
            gamma = self.gamma
            alpha = self.alpha

            first_term = -(p_0 * np.power(self.r(x)/a, phi) / C_11)\
                        * (-np.power(k, 2*n) * (-2*n + S_1) + S_1 * np.power(self.r(x)/a, 2*n))\
                        / ((-1+np.power(k, 2*n)) * (-2*n + S_1) * S_1 * np.power(self.r(x)/a, n))
            second_term = -(p_0 * np.power(self.r(x)/a, phi) / C_11)\
                        * self.F(np.array([1j * gamma]), x) * np.cos(gamma * self.tau)\
                        / self.G(np.array([1j * gamma]))
            series = 0
            for s in range(len(alpha)):
                first_parantheses = S_1 * jv(n, alpha[s] + 0j) - alpha[s] * jv(n-1, alpha[s] + 0j)
                second_parantheses = np.power(S_1, 2) - 2 * S_1 * n + np.power(alpha[s], 2)
                third_parantheses = np.power(S_1 * jv(n, alpha[s] * k + 0j) - alpha[s] * k * jv(n-1, alpha[s] * k + 0j), 2)
                forth_parantheses = S_1 * jv(n, alpha[s] * k + 0j) - alpha[s] * k * jv(n-1, alpha[s] * k + 0j)
                fifth_parantheses = np.power(S_1, 2) - 2 * S_1 * n + np.power(alpha[s] * k, 2)
                sixth_parantheses = np.power(S_1 * jv(n, alpha[s] + 0j) - alpha[s] * jv(n-1, alpha[s] + 0j), 2)

                series += 2 * self.F(np.array([1j * alpha[s]]), x)/(np.power(gamma, 2) - np.power(alpha[s], 2))\
                                        * first_parantheses * forth_parantheses / (second_parantheses * third_parantheses - fifth_parantheses * sixth_parantheses)\
                                        * np.cos(alpha[s] * self.tau)

            third_term = (p_0 * np.power(self.r(x)/a, phi) * np.power(gamma, 2) / C_11) * series

            u = np.zeros((1, x.shape[1]))
            u = first_term + second_term + third_term

            displacement = np.zeros(x.shape)
            mag = np.linalg.norm(x, axis=0)
            for i in range(mag.shape[0]):
                displacement[:, i] = u[i].real * x[:, i] / mag[i]

            if radial==True:
                return u.real
            
            return displacement #u.real # [0]

    radial_disp = []

    u_exact = exact_solution(1, 1, 2, lmbda.value, mu.value, 18/2, gamma_p, num_steps)  

    radial_disp = []
    error = []

    # Stepping through time and solving for displacements and strains
    for i in range(int(num_steps)):
        print(i)
        if i < num_steps:
            p.value = pressure[i]

        # Update the right hand side reusing the initial vector
        with b.localForm() as loc_b:
            loc_b.set(0)
        assemble_vector(b, linear_form)

        # Apply Dirichlet boundary condition to the vector
        apply_lifting(b, [bilinear_form], [[]])
        b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        set_bc(b, [])

        # Solve linear problem
        solver.solve(b, displacement.vector)
        displacement.x.scatter_forward()
        radial_disp.append(np.linalg.norm(displacement.x.array.reshape((-1, 3))[index_of_radius_to_check, :]))

        displacement_exact = fe.fem.Function(V)
        displacement_exact.interpolate(u_exact)
        u_exact.step_in_time()

        error_L2 = MPI.COMM_WORLD.allreduce(
            fe.fem.assemble.assemble_scalar(fe.fem.form(((displacement - displacement_exact))** 2 * ufl.dx))**0.5,
            op=MPI.SUM)


        if mesh.comm.rank == 0:
            error.append(error_L2)

        # Update old fields with new quantities
        update_fields(displacement, old_u, old_velocity, old_acceleration)


    if proc == 0:
        plt.plot(t, [radial_disp[i] for i in range(int(num_steps))])
        plt.show()
        print(error)



if __name__ == "__main__":
    mesh_file = "test.msh"

    main(1.0, 400, 0.3, 18, 0, mesh_file)

Both the plotted radial displacement and the error with respect to the exact solution are different in parallel.
Thank you in advance.

2

I would strongly suggest you try to reduce the complexity of this code, as it is a rather time-consuming effort to read through 350 lines of code to decipher what has gone wrong.

There is also an MPI tutorial for DOLFINx at: Parallel Computations with Dolfinx using MPI — MPI tutorial

3

Thank you @dokken
I read through the tutorial but I couldn’t find out what is causing the problem. I have used apply_lifting , scatter_reverse, scatter_forward, set_bc in the for loop iterating through time at the end of the code as the tutorial.
I tried to reduce the complexity of my code as much as I could.

import dolfinx as fe
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem, assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
import numpy as np
from mpi4py import MPI
import ufl
from petsc4py import PETSc
import matplotlib.pyplot as plt


def main(rho, E, nu, maximum_pressure, minimum_pressure, num_steps=100, element_degree=1):

    ############################################Importing the mesh############################################
    ##########################################################################################################
    mesh = fe.mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([1, 1, 1])],
                         [20, 6, 6], cell_type=fe.mesh.CellType.hexahedron)

    maximum_pressure = fe.default_scalar_type(maximum_pressure)
    minimum_pressure = fe.default_scalar_type(minimum_pressure)
    rho = fe.default_scalar_type(rho)
    mu = fe.fem.Constant(mesh, fe.default_scalar_type(E / (2 * (1 + nu))))
    lmbda = fe.fem.Constant(mesh, fe.default_scalar_type(E * nu / ((1 + nu) * (1 - 2 * nu))))

    ######################################################################################################
    #######################creating function space and boundary condition#################################
    element = ufl.VectorElement("Lagrange", mesh.ufl_cell(), element_degree)
    V = fe.fem.functionspace(mesh, element)

    def clamped_boundary(x):
        return np.isclose(x[0], 0)

    fdim = mesh.topology.dim - 1
    clamped_boundary_facets = fe.mesh.locate_entities_boundary(mesh, fdim, clamped_boundary)
    u_D = np.array([0, 0, 0], dtype=fe.default_scalar_type)
    bc = fe.fem.dirichletbc(u_D, fe.fem.locate_dofs_topological(V, fdim, clamped_boundary_facets), V)

    boundaries = [(1, lambda x: np.isclose(x[0], 1))]

    facet_indices, facet_markers = [], []
    fdim = mesh.topology.dim - 1
    for (marker, locator) in boundaries:
        facets = fe.mesh.locate_entities(mesh, fdim, locator)
        facet_indices.append(facets)
        facet_markers.append(np.full_like(facets, marker))
    facet_indices = np.hstack(facet_indices).astype(np.int32)
    facet_markers = np.hstack(facet_markers).astype(np.int32)
    sorted_facets = np.argsort(facet_indices)
    facet_tag = fe.mesh.meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

    gamma_p = 0.01
    t = np.linspace(0, 2*np.pi / gamma_p, num_steps)
    dt = fe.fem.Constant(mesh, t[1] - t[0])

    pressure = maximum_pressure / 2 * (1 - np.cos(gamma_p * t))
    p = fe.fem.Constant(mesh, pressure[0])   

    ######################################################################################################
    #########################################Formulation##################################################
    u = ufl.TrialFunction(V)           
    v  = ufl.TestFunction(V)          
    old_u  = fe.fem.Function(V)
    old_velocity  = fe.fem.Function(V)
    old_acceleration  = fe.fem.Function(V)

    d = len(u)
    I = ufl.variable(ufl.Identity(d))             # Identity tensor

    metadata = {"quadrature_degree": 2}   
    ds = ufl.Measure('ds', domain=mesh, subdomain_data=facet_tag, metadata=metadata)
    dx = ufl.Measure("dx", domain=mesh, metadata=metadata)

    # Generalized-alpha method parameters
    alpha_m = fe.fem.Constant(mesh, 0.01) 
    alpha_f = fe.fem.Constant(mesh, 0.01) 
    gamma   = fe.fem.Constant(mesh, fe.default_scalar_type(0.5+alpha_f-alpha_m))
    beta    = fe.fem.Constant(mesh, fe.default_scalar_type((gamma+0.5)**2/4.))

    # Update formula for acceleration
    # a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0)
    def update_a(u, u_old, v_old, a_old):
        dt_ = float(dt)
        beta_ = float(beta)
        return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old

    # Update formula for velocity
    # v = dt * ((1-gamma)*a0 + gamma*a) + v0
    def update_v(a, u_old, v_old, a_old):
        dt_ = float(dt)
        gamma_ = float(gamma)
        return v_old + dt_*((1-gamma_)*a_old + gamma_*a)


    def update_fields(u_new, u_old, v_old, a_old):
        u_vec, u0_vec  = u_new.x.array[:], u_old.x.array[:]
        v0_vec, a0_vec = v_old.x.array[:], a_old.x.array[:]

        a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec)
        v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec)

        v_old.x.array[:], a_old.x.array[:] = v_vec, a_vec
        u_old.x.array[:] = u_new.x.array[:]


    def avg(x_old, x_new, alpha):
        return alpha*x_old + (1-alpha)*x_new

    normal = -ufl.FacetNormal(mesh)
    
    def epsilon(u):
        return ufl.sym(ufl.grad(u))

    def S(u):
        return 2.0 * mu * epsilon(u) + lmbda * ufl.nabla_div(u) * I 

    acceleration = update_a(u, old_u, old_velocity, old_acceleration)  
    velocity = update_v(acceleration, old_u, old_velocity, old_acceleration)

    formulation = rho * ufl.inner(avg(old_acceleration, acceleration, alpha_m), v) * dx \
      + ufl.inner(epsilon(v), S(avg(old_u, u, alpha_f))) * dx \
      - ufl.dot(v, p * normal) * ds(1)   

    bilinear_form = fe.fem.form(ufl.lhs(formulation))
    linear_form = fe.fem.form(ufl.rhs(formulation))

    ######################################################################################################
    ###############################################Solving################################################
    A = assemble_matrix(bilinear_form, bcs=[])
    A.assemble()
    b = create_vector(linear_form)

    solver = PETSc.KSP().create(mesh.comm)
    solver.setInitialGuessNonzero(True)
    solver.setOperators(A)
    solver.getPC().setType(PETSc.PC.Type.LU)

    displacement = fe.fem.Function(V)

    class exact_solution():
        def __call__(self, x):
            return x

    u_exact = exact_solution()  
    displacement_exact = fe.fem.Function(V)
    displacement_exact.interpolate(u_exact)

    radial_disp = []
    error = []

    # Stepping through time and solving for displacements
    for i in range(int(num_steps)):
        p.value = pressure[i]

        with b.localForm() as loc_b:
            loc_b.set(0)
        assemble_vector(b, linear_form)

        apply_lifting(b, [bilinear_form], [[]])
        b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        set_bc(b, [])

        solver.solve(b, displacement.vector)
        displacement.x.scatter_forward()
        radial_disp.append(np.linalg.norm(displacement.x.array.reshape((-1, 3))[10, :]))   ###Just a random node

        error_L2 = MPI.COMM_WORLD.allreduce(
            fe.fem.assemble.assemble_scalar(fe.fem.form(((displacement - displacement_exact))** 2 * ufl.dx))**0.5,
            op=MPI.SUM)

        if mesh.comm.rank == 0:
            error.append(error_L2)

        update_fields(displacement, old_u, old_velocity, old_acceleration)


    if MPI.COMM_WORLD.rank == 0:
        plt.plot(t, radial_disp)
        plt.show()
        print(error)


if __name__ == "__main__":

    main(1.0, 400, 0.3, 18, 0)
4

Your L2 error computation is wrong, as \sqrt{a} + \sqrt{b} \neq \sqrt{a+b}.
I.e. you should use:

       error_L2 = MPI.COMM_WORLD.allreduce(
            fe.fem.assemble.assemble_scalar(fe.fem.form(((displacement - displacement_exact))** 2 * ufl.dx)),
            op=MPI.SUM)**0.5

Secondly, for debugging this, I would strongly advice you to check the displacement (visually in Paraview) and check that they look similar in eye-ball norm. That is what I did, and I saw they were the same, but that the error computation gave different results in serial and parallel.

1 Like
5

Thank you really really much for your point,
I have another question about the index of the nodes in parallel, if I want to check for the displacement of one specific node that I know its global index, how can I know which processor with what local index contains the info of that node. Second if I want to save the displacements as petsc binary file like this:

io = PetscBinaryIO.PetscBinaryIO(complexscalars=True)
displacementAsPetscBiIOVec = displacement.vector.array_w.view(PetscBinaryIO.Vec)
io.writeBinaryFile( f"displacement_time_step_{i}.dat", [displacementAsPetscBiIOVec])

Should I append all the displacements calculated by all the processors first?

6

See the post: How to apply a boundary condition to a range of DOFs? - #7 by dokken
on how to map global input node index to a degree of freedom: How to apply a boundary condition to a range of DOFs? - #7 by dokken

For your second question I cannot give many pointers, as I do not use PetscBinaryIO.

1 Like