Create subdomain and integrate over it

[HSadeghein]

I have implemented the lid-driven cavity with a unit square mesh. But I want to extend my mesh to be rectangle and then compute the same thing but only on the bottom half of my mesh. I have used subdomain and updated dx based on my subdomain, but when I call my solve function, it turns a C++ error. malloc(): invalid size (unsorted)

import numpy as np
from numpy import newaxis
from mpi4py import MPI
from dolfinx import mesh
import dolfinx
import pyvista
import dolfinx.plot
from dolfinx.fem import FunctionSpace
from dolfinx import fem
import ufl
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
from tqdm import tqdm
import time

TIME_STEP_LENGTH = 0.01
N_TIME_STEPS = 50
KINEMATIC_VICOSITY = 0.01


def Omega_0(x):
    return np.logical_or(np.isclose(x[1], 1), np.isclose(x[1], 0))


def Omega_1(x):
    return np.logical_and(x[1] <= 1, x[0] <= 1)


def main():
    domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0, 0], [2, 2]], [
                                   40, 40], mesh.CellType.triangle)
    # domain = mesh.create_unit_square(
    #     MPI.COMM_WORLD, 40, 40)

    velocity_function_space = fem.VectorFunctionSpace(domain, ("CG", 2))
    pressure_function_space = fem.FunctionSpace(domain, ("CG", 1))

    u_trial = ufl.TrialFunction(velocity_function_space)
    p_trial = ufl.TrialFunction(pressure_function_space)
    v_test = ufl.TestFunction(velocity_function_space)
    q_test = ufl.TestFunction(pressure_function_space)

    # mf = dolfinx.MeshFunction("size_t", mesh, mesh.topology.dim, 0)
    # top_facets = dolfinx.mesh.locate_entities_boundary(domain, 1, lambda x : np.isclose(x[1], 1))
    # mt = dolfinx.mesh.meshtags(domain, 1, top_facets, 1)
    # dx(1) = ufl.Measure("ds", subdomain_data=mt)
    subdomain = mesh.locate_entities(
        domain, dim=domain.topology.dim, marker=Omega_1)
    subdomain_values = np.full_like(subdomain, 1)
    facet_tags = dolfinx.mesh.meshtags(
        domain, domain.topology.dim, subdomain, subdomain_values)
    dx = ufl.Measure('dx', domain=domain, subdomain_data=facet_tags)
    facets = mesh.locate_entities(domain, dim=(domain.topology.dim - 1),
                                  marker=lambda x: np.logical_or(np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)),
                                                                 np.isclose(x[1], 0.0)))

    facets1 = mesh.locate_entities(domain, dim=(domain.topology.dim - 1),
                                   marker=lambda x: np.isclose(x[1], 1))

    # dx(1) =  ufl.Measure("dx(1)", domain=domain, subdomain_data=facets)
    dofs = fem.locate_dofs_topological(
        V=velocity_function_space, entity_dim=1, entities=facets)

    dofs1 = fem.locate_dofs_topological(
        V=velocity_function_space, entity_dim=1, entities=facets1)
    # u_zero = np.array((0,)*domain.geometry.dim, dtype=ScalarType)
    bc = fem.dirichletbc(ScalarType((0, 0)), dofs=dofs,
                         V=velocity_function_space)
    bc1 = fem.dirichletbc(ScalarType((1, 0)), dofs=dofs1,
                          V=velocity_function_space)

    velocity_boundry_condition = [bc, bc1]

    u_prev = fem.Function(velocity_function_space)
    u_tent = fem.Function(velocity_function_space)
    u_next = fem.Function(velocity_function_space)
    p_next = fem.Function(pressure_function_space)

    momentum_weak_form_residuum = (
        1.0 / TIME_STEP_LENGTH * ufl.inner(u_trial - u_prev, v_test) * dx(1)
        +
        ufl.inner(ufl.grad(u_prev) * u_prev, v_test) * dx(1)
        +
        KINEMATIC_VICOSITY *
        ufl.inner(ufl.grad(u_trial), ufl.grad(v_test)) * dx(1)
    )

    momentum_weak_form_lhs = ufl.lhs(momentum_weak_form_residuum)
    momentum_weak_form_rhs = ufl.rhs(momentum_weak_form_residuum)

    momentum_weak_form_lhs = fem.form(momentum_weak_form_lhs)
    momentum_weak_form_rhs = fem.form(momentum_weak_form_rhs)
    momentum_weak_form_matrix = fem.petsc.assemble_matrix(
        momentum_weak_form_lhs, bcs=velocity_boundry_condition)
    momentum_weak_form_matrix.assemble()
    momentum_weak_form_vector = fem.petsc.create_vector(momentum_weak_form_rhs)

    # Weak form of the pressure
    pressure_possion_weak_lhs = ufl.inner(
        ufl.grad(p_trial), ufl.grad(q_test)) * dx(1)
    pressure_possion_weak_rhs = -1.0 / TIME_STEP_LENGTH * \
        ufl.div(u_tent) * q_test * dx(1)

    pressure_possion_weak_lhs = fem.form(pressure_possion_weak_lhs)
    pressure_possion_weak_rhs = fem.form(pressure_possion_weak_rhs)
    pressure_possion_weak_form_matrix = fem.petsc.assemble_matrix(
        pressure_possion_weak_lhs)
    pressure_possion_weak_form_matrix.assemble()
    pressure_possion_weak_form_vector = fem.petsc.create_vector(
        pressure_possion_weak_rhs)

    # Weak form of the velocity update equation
    velocity_update_weak_form_lhs = ufl.inner(u_trial, v_test) * dx(1)
    velocity_update_weak_form_rhs = ufl.inner(
        u_tent, v_test) * dx(1) - TIME_STEP_LENGTH * ufl.inner(ufl.grad(p_next), v_test) * dx(1)

    velocity_update_weak_form_lhs = fem.form(velocity_update_weak_form_lhs)
    velocity_update_weak_form_rhs = fem.form(velocity_update_weak_form_rhs)
    velocity_update_weak_form_matrix = fem.petsc.assemble_matrix(
        velocity_update_weak_form_lhs, bcs=velocity_boundry_condition)
    velocity_update_weak_form_matrix.assemble()
    velocity_update_weak_form_vector = fem.petsc.create_vector(
        velocity_update_weak_form_rhs)

    # Solver for step 1
    solver1 = PETSc.KSP().create(domain.comm)
    solver1.setOperators(momentum_weak_form_matrix)
    solver1.setType(PETSc.KSP.Type.BCGS)
    pc1 = solver1.getPC()
    pc1.setType(PETSc.PC.Type.HYPRE)
    pc1.setHYPREType("boomeramg")

    # Solver for step 2
    solver2 = PETSc.KSP().create(domain.comm)
    solver2.setOperators(pressure_possion_weak_form_matrix)
    solver2.setType(PETSc.KSP.Type.BCGS)
    pc2 = solver2.getPC()
    pc2.setType(PETSc.PC.Type.HYPRE)
    pc2.setHYPREType("boomeramg")

    # Solver for step 3
    solver3 = PETSc.KSP().create(domain.comm)
    solver3.setOperators(velocity_update_weak_form_matrix)
    solver3.setType(PETSc.KSP.Type.CG)
    pc3 = solver3.getPC()
    pc3.setType(PETSc.PC.Type.SOR)

    OFF_SCREEN = False

    u_plotter = pyvista.Plotter(
        window_size=[1024, 1024], off_screen=OFF_SCREEN)
    u_plotter.open_gif("hassan.gif")

    u_topology, u_cell_types, u_geometry = dolfinx.plot.create_vtk_mesh(domain)
    mesh_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
    u_plotter.add_mesh(mesh_grid, show_edges=True)

    u_topology, u_cell_types, u_geometry = dolfinx.plot.create_vtk_mesh(
        velocity_function_space)
    u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
    values = np.zeros((u_geometry.shape[0], 3), dtype=np.float64)
    values[:, :len(u_next)] = u_next.x.array.real.reshape(
        (u_geometry.shape[0], len(u_next)))
    u_grid["u"] = values
    glyphs = u_grid.glyph(orient="u", factor=0.1)
    u_plotter.add_mesh(glyphs)
    u_plotter.view_xy()
    if not OFF_SCREEN:
        u_plotter.show(interactive_update=True)
    for t in tqdm(range(N_TIME_STEPS)):
        # Step 1: Tentative veolcity step
        with momentum_weak_form_vector.localForm() as loc_1:
            loc_1.set(0)
        fem.petsc.assemble_vector(
            momentum_weak_form_vector, momentum_weak_form_rhs)
        fem.petsc.apply_lifting(momentum_weak_form_vector, [
                                momentum_weak_form_lhs], [velocity_boundry_condition])
        momentum_weak_form_vector.ghostUpdate(
            addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        fem.petsc.set_bc(momentum_weak_form_vector, velocity_boundry_condition)
        solver1.solve(momentum_weak_form_vector, u_tent.vector)
        u_tent.x.scatter_forward()

        # Step 2: Pressure corrrection step
        with pressure_possion_weak_form_vector.localForm() as loc_2:
            loc_2.set(0)
        fem.petsc.assemble_vector(
            pressure_possion_weak_form_vector, pressure_possion_weak_rhs)
        # fem.petsc.apply_lifting(pressure_possion_weak_form_vector, [pressure_possion_weak_lhs])
        # pressure_possion_weak_form_vector.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        # fem.petsc.set_bc(pressure_possion_weak_form_vector)
        solver2.solve(pressure_possion_weak_form_vector, p_next.vector)
        p_next.x.scatter_forward()

        # # Step 3: Velocity correction step
        with velocity_update_weak_form_vector.localForm() as loc_3:
            loc_3.set(0)
        fem.petsc.assemble_vector(
            velocity_update_weak_form_vector, velocity_update_weak_form_rhs)
        velocity_update_weak_form_vector.ghostUpdate(
            addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        solver3.solve(velocity_update_weak_form_vector, u_next.vector)
        u_next.x.scatter_forward()
        # Update variable with solution form this time step
        u_prev.x.array[:] = u_next.x.array[:]
        # p_.x.array[:] = p_.x.array[:]

        # u_plotter.add_mesh(mesh_grid, show_edges=True, color="#222222")
        values = np.zeros((u_geometry.shape[0], 3), dtype=np.float64)
        values[:, :len(u_next)] = u_next.x.array.real.reshape(
            (u_geometry.shape[0], len(u_next)))
        u_grid["u"] = values
        updated_glyphs = u_grid.glyph(orient="u", factor=0.1)

        glyphs.copy_from(updated_glyphs)

        if not OFF_SCREEN:
            u_plotter.add_mesh(glyphs)
        else:
            u_plotter.write_frame()


    u_plotter.close()

if __name__ == "__main__":
    main()

can you please help me with it?

[dokken]

When you define your function space on the whole mesh

@HSadeghein

velocity_function_space = fem.VectorFunctionSpace(domain, ("CG", 2))

you have degrees of freedom for every cell, that means that if you want to deactivate a part of the domain, you need to fix these degrees of freedom in the assembled matrix (otherwise you rows in your matrix with all zeros).

You can create a DirichletBC on all these dofs, to make sure that each row of the matrix with any of these dofs is the identity matrix.

[HSadeghein]

Thank you for your reply!
How should I fix the DOFs on those elements that are unused?

[francesco-ballarin]

You can use for instance multiphenicsx (see https://multiphenics.github.io/ , disclaimer: I am the author of that) or @jpdean branch in dolfinx.

[dokken]

I have shown how to do this in the context of surface integrals at:

jorgensd/dolfinx_mpc/blob/2af8c962f94cf821ea273b568e0336f31cc43513/python/dolfinx_mpc/utils/mpc_utils.py#L47

Which is similar to legacy dolfin:

Obtain velocity normal to boundary

As the FacetNormal as implied by its name, is defined on each Facet, and not for each degree of freedom, you have to approximate the normal at the degrees of freedom (for instance at corners). For your particular case, the problem is that you are trying to interpolate an ComponentTensor, that is an expression using ufl syntax. To work around this, I would suggest using an UserExpression as follows: n = FacetNormal(mesh) # Interpolate mesh values into boundary function u_boundary = Function(D) …

or you could use any of the suggestions by @francesco-ballarin