Skip to content

Commit b64e3b4

Browse files
committed
Fix quadrature rule
1 parent b1ee5c0 commit b64e3b4

File tree

2 files changed

+31
-36
lines changed

2 files changed

+31
-36
lines changed

src/fenicsxconcrete/sensor_definition/stress_sensor.py

Lines changed: 17 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -3,9 +3,8 @@
33
import os
44
from typing import TYPE_CHECKING
55

6-
import numpy as np
76
import dolfinx as df
8-
import ufl
7+
import numpy as np
98

109
if TYPE_CHECKING:
1110
from fenicsxconcrete.finite_element_problem.base_material import MaterialProblem
@@ -43,17 +42,28 @@ def measure(self, problem: MaterialProblem) -> None:
4342
stress_tensor_dim = problem.experiment.mesh.topology.dim
4443
stress_function = project(
4544
stress, # stress fct from problem
46-
df.fem.functionspace(problem.experiment.mesh, (problem.q_fields.plot_space_type[0], problem.q_fields.plot_space_type[1], (stress_tensor_dim,stress_tensor_dim))), # tensor space
45+
df.fem.functionspace(
46+
problem.experiment.mesh,
47+
(
48+
problem.q_fields.plot_space_type[0],
49+
problem.q_fields.plot_space_type[1],
50+
(stress_tensor_dim, stress_tensor_dim),
51+
),
52+
), # tensor space
4753
problem.q_fields.measure,
4854
)
4955
elif mandel_stress is not None:
5056
stress_function = project(
5157
mandel_stress, # stress fct from problem
5258
df.fem.functionspace(
53-
problem.experiment.mesh, (problem.q_fields.plot_space_type[0], problem.q_fields.plot_space_type[1], (problem.mandel_stress_dim,))
54-
), # vector space
55-
ufl.dx, #TODO: check consistency!
56-
#problem.q_fields.measure,
59+
problem.experiment.mesh,
60+
(
61+
problem.q_fields.plot_space_type[0],
62+
problem.q_fields.plot_space_type[1],
63+
(problem.mandel_stress_dim,),
64+
),
65+
),
66+
problem.q_fields.measure,
5767
)
5868
else:
5969
raise Exception("Stress and Mandel stress not defined in problem")

src/fenicsxconcrete/util/quadrature.py

Lines changed: 14 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -40,12 +40,9 @@ def __init__(
4040
self.degree = degree
4141
basix_cell = _ufl_cell_type_to_basix(self.cell_type)
4242
self.points, self.weights = basix.make_quadrature(basix_cell, self.degree, rule=self.type)
43-
self.dx = ufl.dx(
44-
metadata={
45-
"quadrature_rule": self.type,
46-
"quadrature_degree": self.degree,
47-
}
48-
)
43+
44+
metadata = {"quadrature_degree": self.degree, "quadrature_scheme": self.type.name}
45+
self.dx = ufl.dx(metadata=metadata)
4946

5047
def create_quadrature_space(self, mesh: df.mesh.Mesh) -> df.fem.functionspace:
5148
"""
@@ -56,13 +53,9 @@ def create_quadrature_space(self, mesh: df.mesh.Mesh) -> df.fem.functionspace:
5653
A scalar quadrature `FunctionSpace` on `mesh`.
5754
"""
5855
assert mesh.ufl_cell() == self.cell_type
59-
Qe = basix.ufl.quadrature_element(mesh.topology.cell_name(), value_shape=(), degree=self.degree)
60-
# Qe = ufl.FiniteElement(
61-
# "Quadrature",
62-
# self.cell_type,
63-
# self.degree,
64-
# quad_scheme=self.type.name,
65-
# )
56+
Qe = basix.ufl.quadrature_element(
57+
mesh.topology.cell_name(), value_shape=(), degree=self.degree, scheme=self.type.name
58+
)
6659

6760
return df.fem.functionspace(mesh, Qe)
6861

@@ -76,14 +69,10 @@ def create_quadrature_vector_space(self, mesh: df.mesh.Mesh, dim: int) -> df.fem
7669
A vector valued quadrature `FunctionSpace` on `mesh`.
7770
"""
7871
assert mesh.ufl_cell() == self.cell_type
79-
# Qe = ufl.VectorElement(
80-
# "Quadrature",
81-
# self.cell_type,
82-
# self.degree,
83-
# quad_scheme=self.type.name,
84-
# dim=dim,
85-
# )
86-
Qe = basix.ufl.quadrature_element(mesh.topology.cell_name(), value_shape=(dim,), degree=self.degree)
72+
73+
Qe = basix.ufl.quadrature_element(
74+
mesh.topology.cell_name(), value_shape=(dim,), degree=self.degree, scheme=self.type.name
75+
)
8776

8877
return df.fem.functionspace(mesh, Qe)
8978

@@ -97,14 +86,10 @@ def create_quadrature_tensor_space(self, mesh: df.mesh.Mesh, shape: tuple[int, i
9786
A tensor valued quadrature `FunctionSpace` on `mesh`.
9887
"""
9988
assert mesh.ufl_cell() == self.cell_type
100-
# Qe = ufl.TensorElement(
101-
# "Quadrature",
102-
# self.cell_type,
103-
# self.degree,
104-
# quad_scheme=self.type.name,
105-
# shape=shape,
106-
# )
107-
Qe = basix.ufl.quadrature_element(mesh.topology.cell_name(), value_shape=shape, degree=self.degree)
89+
90+
Qe = basix.ufl.quadrature_element(
91+
mesh.topology.cell_name(), value_shape=shape, degree=self.degree, scheme=self.type.name
92+
)
10893

10994
return df.fem.functionspace(mesh, Qe)
11095

0 commit comments

Comments
 (0)