Hi all,
What is the equivalent code if Q is a MixedSpace and I wish to assign the values (like below) only to a function defined on Q.sub(1) subspace.
Code taken from here.
Q = FunctionSpace(mesh, ("DG", 0))
kappa = Function(Q)
bottom_cells = ct.find(bottom_marker)
kappa.x.array[bottom_cells] = np.full_like(bottom_cells, 1, dtype=ScalarType)
top_cells = ct.find(top_marker)
kappa.x.array[top_cells] = np.full_like(top_cells, 0.1, dtype=ScalarType)
Thanks for your time.
Hi @dokken ,
Perhaps, I should have been more clear with my query. What I need is to be able to assign values to a function using a cell tag like “ct.find(12)”, where the integer 12 is an identifier/marker of a region in the mesh prepared using GMSH’s physical surface utility.
As of now this gives results in an error:
kappa = fem.Function(V=MixedSpace)
W, Wmap = MixedSpace.sub(1).collapse()
bottom_cells = ct.find(14)
kappa.x.array[Wmap] = np.full_like(bottom_cells, 1, dtype=PETSc.ScalarType)
top_cells = ct.find(13)
kappa.x.array[Wmap] = np.full_like(top_cells, 0.1, dtype=PETSc.ScalarType)
Please provide a reproducible code, so that I can see the actual error message and the cause of it.
You are missing a step here, as you need to actually assign the data (ie. 1 and 0.1) to a function in W, and then assign that array to kappa.
If you dont want to do that, you need to use the correct subset of Wmap, i.e. Wmap[bottom_cells] and Wmap[top_cells] in your assignment.
Hi @dokken ,
The mesh file is located here in zenodo.
And here is my attempt:
import numpy as np
from petsc4py import PETSc
import ufl
from ufl import (MixedElement)
from dolfinx import (fem, io)
from dolfinx.io.gmshio import read_from_msh
from mpi4py import MPI
"""
External mesh loading
"""
domain, ct, ft = read_from_msh(
filename="mesh.msh",
comm=MPI.COMM_WORLD,
rank=0, gdim=2
)
gdim = domain.topology.dim
fdim = gdim - 1
n = ufl.FacetNormal(
domain=domain
)
surf_upside = 14
surf_dnside = 13
line_noslip = 12
"""
Spaces and Functions
"""
degree = 2
DG_elem = ufl.VectorElement(
family="Discontinuous Lagrange",
cell=domain.ufl_cell(),
degree=degree
)
CG_elem = ufl.FiniteElement(
family="Lagrange",
cell=domain.ufl_cell(),
degree=degree - 1
)
MixedSpace = fem.FunctionSpace(
mesh=domain,
element=MixedElement([DG_elem, CG_elem])
)
v, q = ufl.TestFunctions(
function_space=MixedSpace
)
state = fem.Function(
V=MixedSpace,
name="Current state"
)
"""
Load and assign the initial condition
"""
# def init_h(x):
# down_idx = x[1] > 0.4
# values = np.full((x.shape[1],), 0.3)
# values[down_idx] = 0.01
# return values
#
#
# state.sub(1).interpolate(lambda x: init_h(x))
space1, map1 = state.sub(1).collapse()
bottom_cells = ct.find(13)
space1[map1[bottom_cells]] = np.full_like(map1[bottom_cells], 1, dtype=PETSc.ScalarType)
top_cells = ct.find(14)
space1[map1[top_cells]] = np.full_like(map1[top_cells], 0.1, dtype=PETSc.ScalarType)
with io.XDMFFile(comm=MPI.COMM_WORLD,
filename="state.xdmf",
file_mode="w",
encoding=io.XDMFFile.Encoding.HDF5) as xdmf:
xdmf.write_mesh(mesh=domain)
xdmf.write_function(u=state.sub(1).collapse())
In order to assign an initial condition to a variable that lies in the second subspace of the mixed space, I have specified two physical surfaces in GMSH. Orange part denotes cells with tag 13. Green part denotes cells with tag 14. Therefore, in the code, I wish to load the mesh, identify the regions with cell tags and assign different values to them. I wish to achieve the result like the commented code (when uncommented) but with cell tags and not with geometrical condition.
Thank you for your patience!
This link does not work.
Also, you are saying that you want to set this to the second space in the mixed space, but that is a CG-1 function, not a DG-0 function as the original question related to.
Here is a MWE that produces what you want.
Make sure you spend some time considering the differences between your code and this one:
import numpy as np
from petsc4py import PETSc
import ufl
from dolfinx import (fem, io, mesh)
from mpi4py import MPI
domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
one_side = mesh.locate_entities(
domain, domain.topology.dim, lambda x: np.logical_and(x[0] < 0.5, x[1] < 0.75))
all_cells_on_procs = domain.topology.index_map(
domain.topology.dim).size_local + domain.topology.index_map(domain.topology.dim).num_ghosts
cells = np.arange(all_cells_on_procs, dtype=np.int32)
surf_upside = 14
surf_dnside = 13
markers = np.full_like(cells, surf_upside, dtype=PETSc.ScalarType)
markers[one_side] = surf_dnside
ct = mesh.meshtags(domain, domain.topology.dim, cells, markers)
degree = 2
DG_elem = ufl.VectorElement(
family="Discontinuous Lagrange",
cell=domain.ufl_cell(),
degree=degree
)
marker_el = ufl.FiniteElement(
family="Discontinuous Lagrange",
cell=domain.ufl_cell(), degree=0)
MixedSpace = fem.FunctionSpace(
mesh=domain,
element=ufl.MixedElement([DG_elem, marker_el])
)
state = fem.Function(MixedSpace, name="Current_state")
space1, map1 = MixedSpace.sub(1).collapse()
map1 = np.asarray(map1, dtype=np.int32)
bottom_cells = ct.find(surf_dnside)
state.x.array[map1[bottom_cells]] = np.full_like(
map1[bottom_cells], 1., dtype=PETSc.ScalarType)
top_cells = ct.find(surf_upside)
state.x.array[map1[top_cells]] = np.full_like(
map1[top_cells], 0.1, dtype=PETSc.ScalarType)
with io.XDMFFile(comm=MPI.COMM_WORLD,
filename="state.xdmf",
file_mode="w",
encoding=io.XDMFFile.Encoding.HDF5) as xdmf:
xdmf.write_mesh(mesh=domain)
xdmf.write_function(u=state.sub(1).collapse())
Dear @dokken ,
Your help is invaluable.
The code solved my problem.
Sorry, the page had not been published previously. The mesh can be found here now.
Thanks again!
I think I understand the difference but for degree=2 and “Lagrange” for the second space, why does the code give error. For higher degrees it works. Am I missing something!
import numpy as np
from petsc4py import PETSc
import ufl
from dolfinx import (fem, io, mesh)
from mpi4py import MPI
from dolfinx.io.gmshio import read_from_msh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
one_side = mesh.locate_entities(
domain, domain.topology.dim, lambda x: np.logical_and(x[0] < 0.5, x[1] < 0.75))
all_cells_on_procs = domain.topology.index_map(
domain.topology.dim).size_local + domain.topology.index_map(domain.topology.dim).num_ghosts
cells = np.arange(all_cells_on_procs, dtype=np.int32)
surf_upside = 14
surf_dnside = 13
markers = np.full_like(cells, surf_upside, dtype=PETSc.ScalarType)
markers[one_side] = surf_dnside
ct = mesh.meshtags(domain, domain.topology.dim, cells, markers)
degree = 2
DG_elem = ufl.VectorElement(
family="Discontinuous Lagrange",
cell=domain.ufl_cell(),
degree=degree
)
CG_elem = ufl.FiniteElement(
family="Lagrange",
cell=domain.ufl_cell(),
degree=degree-1
)
MixedSpace = fem.FunctionSpace(
mesh=domain,
element=ufl.MixedElement([DG_elem, CG_elem])
)
state = fem.Function(MixedSpace, name="Current_state")
space1, map1 = MixedSpace.sub(1).collapse()
map1 = np.asarray(map1, dtype=np.int32)
bottom_cells = ct.find(surf_dnside)
state.x.array[map1[bottom_cells]] = np.full_like(
map1[bottom_cells], 1., dtype=PETSc.ScalarType)
top_cells = ct.find(surf_upside)
state.x.array[map1[top_cells]] = np.full_like(
map1[top_cells], 0.1, dtype=PETSc.ScalarType)
with io.XDMFFile(comm=MPI.COMM_WORLD,
filename="state.xdmf",
file_mode="w",
encoding=io.XDMFFile.Encoding.HDF5) as xdmf:
xdmf.write_mesh(mesh=domain)
xdmf.write_function(u=state.sub(1).collapse())
You are trying to assign data based on the cell indices.
Then you need a function-space that is one-to-one with the cell indices. Only DG-0 is this, as it is a constant per element. Using any other function space to get data that is defined as a constant per cell does not make sense.
You are right!
I need to be careful while assigning the values with this method.
Thanks for your time.