Possible to assemble load vector only for some nodes?

Hi :slight_smile:

Consider a Poisson equation u’’ = f with its assembled load vector (the right-hand side) L:

domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) 
Vh = fem.FunctionSpace(domain, ("CG", 1))
v = ufl.TestFunction(Vh)
u = ufl.TrialFunction(Vh)
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * v
L = fem.assemble_vector(fem.form(f * v * ufl.dx))

This vector L has N lines, where N is the number of nodes in the mesh. I want to implement DEIM using dolfinx. I won’t go into details of DEIM, but it requires computing only some lines of the load vector L (to speed up the computation) rather than assembling for all N nodes. Is there any way I can assemble only for some nodes?

Of course, I can construct L and then select some of its rows, but it kills the performance gain the DEIM method is about. I really appreciate your help/ideas.

If you can link the nodes to certain cell indices, you can restrict the integration domain by using a MeshTag as input to the dx-Measure, as shown in: Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial (for facet integrals, but the same principle works for cell integrals), see for instance: dolfinx/test_assemble_domains.py at 2be366d20632377e0688df585fcd7e201d5dc86c · FEniCS/dolfinx · GitHub or for instance Add PML demo - GSoC 2022 by mikics · Pull Request #2276 · FEniCS/dolfinx · GitHub

1 Like

Thank you! That was so helpful to understand the dolfinx possibilities. I tried to build upon this one that you suggested but I have " incompatible constructor arguments" error for fem.form. Do you have any idea?

import ufl
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem
import pytest
domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

Vh = fem.FunctionSpace(domain, ("CG", 1))
v = ufl.TestFunction(Vh)
u = ufl.TrialFunction(Vh)

cell_map = domain.topology.index_map(domain.topology.dim)
num_cells = cell_map.size_local
indices = np.arange(0, num_cells)
values = np.full(indices.shape, 3, dtype=np.intc)
values[0] = 1
values[1] = 2

marker = pytest.mark.meshtags_factory(domain, domain.topology.dim, indices, values)
dz = ufl.Measure('dx', domain=domain, subdomain_data=marker)
a = fem.form(v * dz(3))

Try:

import ufl
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem
import pytest
domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

Vh = fem.FunctionSpace(domain, ("CG", 1))
v = ufl.TestFunction(Vh)
u = ufl.TrialFunction(Vh)

cell_map = domain.topology.index_map(domain.topology.dim)
num_cells = cell_map.size_local
indices = np.arange(0, num_cells)
values = np.full(indices.shape, 3, dtype=np.intc)
values[0] = 1
values[1] = 2

marker = mesh.meshtags(domain, domain.topology.dim, indices, values)
dz = ufl.Measure('dx', domain=domain, subdomain_data=marker)
a = fem.form(v * dz(3))
1 Like

Great! Thanks a lot. It is working well now. My only surprise is that the assembling even got slower using this marker (and dz(3) rather than dz), and I am wondering why.

Please make a minimal example where you show timings of this claim.

I take back what I said, I was not precise. The issue is not about assembling time, I used it inside a nonlinear solver and the timing I have is not directly related to this dz measure, but lots of other stuff in convergence for example. Anyway, thank you, my code is at least working :slight_smile:

Hi, I have a following question here: how to define L and solver in this case? I followed your steps and it reports the error by petsc-solver: object is in wrong state/Matrix is missing diagonal entry

If you want to solve a problem on only a part of your mesh, you can either:

  1. add Dirichlet BCs on all remaining dofs (say, set then to 0).
  2. use a sub mesh and solve the problem with a function space defined on the sub mesh.

Hi dokken, thank you for your prompt answer.
To 2: Do you have a example of how to define a sub mesh?

See:

Thanks. I’ll look through it. I would like to ask one more question if it’s relevant: I now have a 3D beam mesh and I want to give it a torsion. What may be the most efficient way to do this in fenicsx?

See for instance: dolfinx/main.cpp at main · FEniCS/dolfinx · GitHub
which is similar to: https://fenicsproject.org/olddocs/dolfin/latest/python/demos/hyperelasticity/demo_hyperelasticity.py.html in legacy FEniCS.

Much Thanks. I think it will also help other people. C++ is for now too hard for me to understand. But firslty I will try to implement the old code from Fenics.