Questions about real elements in scifem

General
1

Hi,

I need to use real elements from scifem libary and I would appreciate if someone could answer a few questions I have. The relevant solved example is this:

  1. I noticed that in defining the weak form, simply replacing ufl.dx with dx = ufl.Measure('dx', domain=mesh, subdomain_data=cell_tags) leads to a long cryptic error in the bilinear form compilation (and perhaps in linear too). Is it really prevented to use ufl.Measure('dx') when formulating the problem this way or there is a correct way to do it?

  2. How the code should be modified when the problem is non-linear? Any examples I could study?

Thanks a lot!

2

Without an example of your code and the error, it’s very hard to help you.

The real-function space is a function space as any other space in dolfinx, and you could for instance use the blockedNewtonSolver of scifem: Nonlinear elasticity with blocked Newton solver — scifem
or
scifem/tests/test_blocked_newton_solver.py at main · scientificcomputing/scifem · GitHub

1 Like
3

I thought I had closely adapted the example model so an MWE was redundant. But it is embarrassing to acknowledge how I still terribly underestimate the significance of this step :slight_smile:

While turning an MWE, I realized that I can indeed ufl.Measure('dx'). At least in my specific case, I had a sum of multiple terms in one expression and I was using ufl.dx with some and ufl.Measure(‘dx’) with others. Uniformly using ufl.Measure('dx') with all terms in that line has solved the problem and I will happily take it for now

Thanks a lot for your help as always!

4

I think this is resolved on the main branch as the following PR was merged:

which linked to something very similar:

1 Like
5

I am taking your word that real function space is like any other and trying to solve it how I would ordinarily solve a mixed function space problem. The motivation is to use ufl.lhs and ufl.rhs to automatically extract the bilinear and linear forms. I am expecting to have an array of real elements which will make it terrible to manually construct the block matrix.

However, I am getting an error with ufl.lhs and ufl.rhs operations that I do not understand. I would really appreciate if you could fix the adaptation below of the example model

Thanks a lot!

EDIT: I am using dolfinx 0.9 through conda

from packaging.version import Version
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.cpp.la.petsc import scatter_local_vectors, get_local_vectors
import dolfinx.fem.petsc

import numpy as np
from scifem import create_real_functionspace, assemble_scalar
import ufl

M = 20
mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, M, M, dolfinx.mesh.CellType.triangle, dtype=np.float64
)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))


def u_exact(x):
    return 0.3 * x[1] ** 2 + ufl.sin(2 * ufl.pi * x[0])


x = ufl.SpatialCoordinate(mesh)
n = ufl.FacetNormal(mesh)
g = ufl.dot(ufl.grad(u_exact(x)), n)
f = -ufl.div(ufl.grad(u_exact(x)))
h = assemble_scalar(u_exact(x) * ufl.dx)

# We then create the Lagrange multiplier space

R = create_real_functionspace(mesh)

# Next, we can create a mixed-function space for our problem

if dolfinx.__version__ == "0.8.0":
    u = ufl.TrialFunction(V)
    lmbda = ufl.TrialFunction(R)
    du = ufl.TestFunction(V)
    dl = ufl.TestFunction(R)
elif Version(dolfinx.__version__) >= Version("0.9.0.0"):
    W = ufl.MixedFunctionSpace(V, R)
    u, lmbda = ufl.TrialFunctions(W)
    du, dl = ufl.TestFunctions(W)
else:
    raise RuntimeError("Unsupported version of dolfinx")

# We can now define the variational problem

zero = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0))

dx = ufl.Measure('dx', mesh)

F = ufl.inner(ufl.grad(u), ufl.grad(du)) * ufl.dx
F += ufl.inner(lmbda, du) * ufl.dx
F += ufl.inner(u, dl) * ufl.dx
F += -ufl.inner(f, du) * ufl.dx + ufl.inner(g, du) * ufl.ds
F += -ufl.inner(zero, dl) * ufl.dx

a = ufl.lhs(F)
L = ufl.rhs(F)

from dolfinx.fem.petsc import LinearProblem

petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
problem = LinearProblem(a, L, bcs=[], petsc_options=petsc_options)
uh = problem.solve()