Skip to content

Commit

Permalink
Possible fix for taking scalar derivatives of vector coefficients usi…
Browse files Browse the repository at this point in the history
…ng coefficient derivative maps (#124)

* A possible fix for taking the derivative of vector coefficients w.r.t. a scalar.

Also adding a test.

* A slightly neater version that handles scalar derivatives of vectors.

---------

Co-authored-by: Matthew Scroggs <[email protected]>
  • Loading branch information
cianwilson and mscroggs authored Sep 12, 2023
1 parent 551f673 commit 37ae353
Show file tree
Hide file tree
Showing 2 changed files with 27 additions and 3 deletions.
24 changes: 24 additions & 0 deletions test/test_derivative.py
Original file line number Diff line number Diff line change
Expand Up @@ -483,6 +483,30 @@ def test_coefficient_derivatives(self):
self.assertEqual(replace(actual, fd.function_replace_map), expected)


def test_vector_coefficient_scalar_derivatives(self):
V = FiniteElement("Lagrange", triangle, 1)
VV = VectorElement("Lagrange", triangle, 1)

dv = TestFunction(V)

df = Coefficient(VV, count=0)
g = Coefficient(VV, count=1)
f = Coefficient(VV, count=2)
u = Coefficient(V, count=3)
cd = {f: df}

integrand = inner(f, g)

i0, i1, i2, i3, i4 = [Index(count=c) for c in range(5)]
expected = as_tensor(df[i1]*dv, (i1,))[i0]*g[i0]

F = integrand*dx
J = derivative(F, u, dv, cd)
fd = compute_form_data(J)
actual = fd.preprocessed_form.integrals()[0].integrand()
assert (actual*dx).signature() == (expected*dx).signature()


def test_vector_coefficient_derivatives(self):
V = VectorElement("Lagrange", triangle, 1)
VV = TensorElement("Lagrange", triangle, 1)
Expand Down
6 changes: 3 additions & 3 deletions ufl/algorithms/apply_derivatives.py
Original file line number Diff line number Diff line change
Expand Up @@ -835,9 +835,9 @@ def coefficient(self, o):
dosum = Zero(o.ufl_shape)
for do, v in zip(dos, self._v):
so, oi = as_scalar(do)
rv = len(v.ufl_shape)
oi1 = oi[:-rv]
oi2 = oi[-rv:]
rv = len(oi) - len(v.ufl_shape)
oi1 = oi[:rv]
oi2 = oi[rv:]
prod = so * v[oi2]
if oi1:
dosum += as_tensor(prod, oi1)
Expand Down

0 comments on commit 37ae353

Please sign in to comment.