# Derive weak form from variational principle
**URL:** https://fenicsproject.discourse.group/t/derive-weak-form-from-variational-principle/8894
**Category:** variational formulation
**Created:** [August 2, 2022, 12:16pm UTC](https://fenicsproject.discourse.group/t/derive-weak-form-from-variational-principle/8894 "2022-08-02T12:16:03Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@pamgur](https://fenicsproject.discourse.group/u/pamgur)
#### Post date: [August 2, 2022, 12:16pm UTC](https://fenicsproject.discourse.group/t/derive-weak-form-from-variational-principle/8894/1 "2022-08-02T12:16:03Z")
Hi all,
I’m trying to solve a set of coupled non-linear PDEs in FEniCSx, but I’m having trouble understanding how to write the variational problem. For now I’m assuming Dirichlet BCs, although later I will need to include rather complicated Robin BCs. Moreover, the problem is such that spherical coordinates and basis are the natural set to work on.
First of all, we have an energy functional F that depends on the vector order parameter \mathbf{m}\equiv(m\_r,m\_{\theta},m\_{\phi}) and on its derivatives, which for simplicity I note as \partial\mathbf{m}:
F[\mathbf{m},\partial\mathbf{m}]= \int\_{\Omega} \mathrm{d}V f\_B(\mathbf{m},\partial\mathbf{m}) = \int\_{\Omega} \mathrm{d}r \mathrm{d}\theta \mathrm{d}\phi\, J(r, \theta) f\_B(\mathbf{m},\partial\mathbf{m})
with J(r, \theta)=r^2\sin\theta. Without specifying f\_B, I can impose \delta F =0 to minimise the energy, and thus derive the Euler-Lagrange equations for m\_i,\ i=r,\theta,\phi:
J(r, \theta)\left[\dfrac{\partial f\_B}{\partial m\_i}-\nabla\cdot \dfrac{\partial f\_B}{\partial (\partial m\_i)}\right]=0
where \nabla\cdot is the divergence in spherical coordinates and
\dfrac{\partial f\_B}{\partial (\partial m\_i)}\equiv \left( \dfrac{\partial f\_B}{\partial (\partial\_r m\_i)}, \dfrac{\partial f\_B}{\partial (\partial\_{\theta} m\_i)},\dfrac {\partial f\_B}{\partial (\partial\_{\phi} m\_i)}\right)
is a vector containing the derivatives of f\_B with respect to the derivatives of m\_i.
\*\*\*\*\*\ ***First question:** can I derive the weak formulation directly from this general equation or (even better) from the energy functional, and then replace with f\_B? Or is it necessary to specify f\_B from the beginning?
In this case,
f\_B=\dfrac{a}{2}m^2+\dfrac{c}{4}m^4 + d \left\{\left(\nabla\cdot\mathbf m\right)^2 + \left[\mathbf m\cdot \left(\nabla\times\mathbf m\right)\right]^2 + \left[\mathbf m\times \left(\nabla\times\mathbf m\right)\right]^2\right\} \,.
Despite assuming rotational symmetry (so nothing depends on \phi) and taking m\_{\phi}=0, this gives rise to very complicated equations for the remaining components of \mathbf{m}\equiv(m\_r,m\_{\theta},0):
(As an image from Mathematica so I don’t make any mistakes while copying them.)
\*\*\*\*\*\ ***Second question:** Any suggestions on how to derive the weak form efficiently and then implement it in FEniCSx?
Thanks in advance!
---
### Author:  [@kamensky](https://fenicsproject.discourse.group/u/kamensky)
#### Post date: [August 2, 2022, 4:58pm UTC](https://fenicsproject.discourse.group/t/derive-weak-form-from-variational-principle/8894/2 "2022-08-02T16:58:46Z")
The canonical approach in FEniCS would be to specify F in UFL, then use the `derivative()` function to automatically take its Gateaux derivative with respect to \mathbf{m} in the direction of a test function \delta\mathbf{m} to get \delta F, which would then be the residual of the weak problem. This is done in the Legacy FEniCS hyperelasticity [demo](https://fenicsproject.org/olddocs/dolfin/latest/python/demos/hyperelasticity/demo_hyperelasticity.py.html), which should be easy to adapt to FEniCSx. (The current FEniCSx hyperelasticity [demo](https://jorgensd.github.io/dolfinx-tutorial/chapter2/hyperelasticity.html) formulates the problem in a different but equivalent way, by taking a partial derivative of the energy density with `diff()`, to get the first Piola–Kirchhoff stress; this is slightly more work, but can be clearer when combining with time integration and non-conservative forcing. That would be analogous to starting from your Euler–Lagrange form, integrating by parts, and only automating the partials of f\_B with respect to \mathbf{m} and \nabla\mathbf{m}.)
To handle the change of coordinates, it would probably be cleanest to define auxiliary Python functions for the polar divergence and curl, interpreting UFL’s `SpatialCoordinate` to be (r,\theta,\phi), where UFL `grad()` would be a vector of partial derivatives with respect to polar coordinates (and not a true coordinate-independent gradient) and a `dx`-type `Measure` would similarly be dr\,d\theta\,d\phi (instead of a physical volume measure).
---
### Author:  [@pamgur](https://fenicsproject.discourse.group/u/pamgur)
#### Post date: [August 3, 2022, 10:20am UTC](https://fenicsproject.discourse.group/t/derive-weak-form-from-variational-principle/8894/3 "2022-08-03T10:20:31Z")
Hi kamensky, thanks for your reply. I’ve started by defining functions for the spherical curl, divergence and gradient, as you suggested, but I don’t know how to return a UFL type instead of a list or numpy array. As a consequence, I can’t compute dot products between a vector function and its curl, for example. Here’s a mwe:
```auto
import ufl
from dolfinx import mesh, fem
from mpi4py import MPI
dim = 3 # dimension of the OP
domain = mesh.create_unit_cube(MPI.COMM_WORLD,10,10,10)
x = ufl.SpatialCoordinate(domain)
element = ufl.VectorElement("CG", domain.ufl_cell(), 1, dim)
V = fem.FunctionSpace(domain, element) # space of vector functions with dim components
m = ufl.TrialFunction(V) # dim-dimensional vector OP
def curl_sph(A, x):
xr, xq, xf = x; jr, jq, jf = 0, 1, 2
Ar, Aq, Af = A
curl_r = (ufl.Dx(Af*ufl.sin(xq),jq)-ufl.Dx(Aq,jf))/(xr*ufl.sin(xq))
curl_q = (ufl.Dx(Ar,jf)/ufl.sin(xq)-ufl.Dx(xr*Af,jr))/xr
curl_f = (ufl.Dx(xr*Aq,jr)-ufl.Dx(Ar,jq))/xr
return [curl_r, curl_q, curl_f]
#ufl.dot([1,0,0],[0,1,0])
ufl.dot(m,curl_sph(m,x))
```
```auto
---------------------------------------------------------------------------
UFLValueError Traceback (most recent call last)
Input In [1], in ()
20 return [curl_r, curl_q, curl_f]
22 #ufl.dot([1,0,0],[0,1,0])
---> 23 ufl.dot(m,curl_sph(m,x))
File /usr/local/lib/python3.10/dist-packages/ufl/operators.py:173, in dot(a, b)
171 "UFL operator: Take the dot product of *a* and *b*. This won't take the complex conjugate of the second argument."
172 a = as_ufl(a)
--> 173 b = as_ufl(b)
174 if a.ufl_shape == () and b.ufl_shape == ():
175 return a * b
File /usr/local/lib/python3.10/dist-packages/ufl/constantvalue.py:437, in as_ufl(expression)
435 return IntValue(expression)
436 else:
--> 437 raise UFLValueError("Invalid type conversion: %s can not be converted"
438 " to any UFL type." % str(expression))
UFLValueError: Invalid type conversion: [Division(Sum(Indexed(Grad(Product(Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None), MultiIndex((FixedIndex(2),))), Sin(Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(1),)))))), MultiIndex((FixedIndex(1),))), Product(IntValue(-1), Indexed(Grad(Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None), MultiIndex((FixedIndex(1),)))), MultiIndex((FixedIndex(2),))))), Product(Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(0),))), Sin(Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(1),)))))), Division(Sum(Product(IntValue(-1), Indexed(Grad(Product(Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(0),))), Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None), MultiIndex((FixedIndex(2),))))), MultiIndex((FixedIndex(0),)))), Division(Indexed(Grad(Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None), MultiIndex((FixedIndex(0),)))), MultiIndex((FixedIndex(2),))), Sin(Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(1),)))))), Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(0),)))), Division(Sum(Indexed(Grad(Product(Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(0),))), Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None), MultiIndex((FixedIndex(1),))))), MultiIndex((FixedIndex(0),))), Product(IntValue(-1), Indexed(Grad(Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None), MultiIndex((FixedIndex(0),)))), MultiIndex((FixedIndex(1),))))), Indexed(SpatialCoordinate(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0)), MultiIndex((FixedIndex(0),))))] can not be converted to any UFL type.
```
In fact, this extremely minimal example raises the same error:
```auto
import ufl
ufl.dot([1,0,0],[0,1,0])
```
```auto
---------------------------------------------------------------------------
UFLValueError Traceback (most recent call last)
Input In [1], in | ()
1 import ufl
----> 3 ufl.dot([1,0,0],[0,1,0])
File /usr/local/lib/python3.10/dist-packages/ufl/operators.py:172, in dot(a, b)
170 def dot(a, b):
171 "UFL operator: Take the dot product of *a* and *b*. This won't take the complex conjugate of the second argument."
--> 172 a = as_ufl(a)
173 b = as_ufl(b)
174 if a.ufl_shape == () and b.ufl_shape == ():
File /usr/local/lib/python3.10/dist-packages/ufl/constantvalue.py:437, in as_ufl(expression)
435 return IntValue(expression)
436 else:
--> 437 raise UFLValueError("Invalid type conversion: %s can not be converted"
438 " to any UFL type." % str(expression))
UFLValueError: Invalid type conversion: [1, 0, 0] can not be converted to any UFL type.
```
How can I return something compatible with ufl.dot?
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [August 3, 2022, 10:35am UTC](https://fenicsproject.discourse.group/t/derive-weak-form-from-variational-principle/8894/4 "2022-08-03T10:35:48Z")
> [@pamgur](#):
>
> ` return [curl_r, curl_q, curl_f]`
Use `ufl.as_vector([curl_r, curl_q, curl_f])`.
---
### Author:  [@pamgur](https://fenicsproject.discourse.group/u/pamgur)
#### Post date: [August 4, 2022, 9:18am UTC](https://fenicsproject.discourse.group/t/derive-weak-form-from-variational-principle/8894/5 "2022-08-04T09:18:40Z")
Hi again. So thanks to @dokken I’ve been able to write the energy, but now I can’t take the Gateaux derivative as suggested by @kamensky. I’ve followed the steps in [Hyperelasticity](https://fenicsproject.org/olddocs/dolfin/latest/python/demos/hyperelasticity/demo_hyperelasticity.py.html). Here’s a MWE:
```auto
import ufl
from dolfinx import mesh, fem
from mpi4py import MPI
dim = 3 # dimension of the OP
domain = mesh.create_unit_cube(MPI.COMM_WORLD,10,10,10)
x = ufl.SpatialCoordinate(domain)
element = ufl.VectorElement("CG", domain.ufl_cell(), 1, dim)
V = fem.FunctionSpace(domain, element) # space of vector functions with dim components
def Jacob(x):
"""Jacobian of the transformation from Cartesian to spherical coordinates."""
xr, xq, xf = x
return xr**2*ufl.sin(xq)
m = ufl.TrialFunction(V) # dim-dimensional vector OP
phi = ufl.TestFunction(V) # test function
F = ufl.dot(m,m)*Jacob(x)*ufl.dx
dF = ufl.derivative(F, m, phi) # variation of the energy
```
which returns the error:
```auto
ERROR:UFL:Invalid coefficient type for
---------------------------------------------------------------------------
UFLException Traceback (most recent call last)
Input In [2], in ()
20 phi = ufl.TestFunction(V) # test function
22 F = ufl.dot(m,m)*Jacob(x)*ufl.dx
---> 23 dF = ufl.derivative(F, m, phi)
File /usr/local/lib/python3.10/dist-packages/ufl/formoperators.py:272, in derivative(form, coefficient, argument, coefficient_derivatives)
248 def derivative(form, coefficient, argument=None, coefficient_derivatives=None):
249 """UFL form operator:
250 Compute the Gateaux derivative of *form* w.r.t. *coefficient* in direction
251 of *argument*.
(...)
269 ``Coefficient`` instances to their derivatives w.r.t. *coefficient*.
270 """
--> 272 coefficients, arguments = _handle_derivative_arguments(form, coefficient,
273 argument)
275 if coefficient_derivatives is None:
276 coefficient_derivatives = ExprMapping()
File /usr/local/lib/python3.10/dist-packages/ufl/formoperators.py:222, in _handle_derivative_arguments(form, coefficient, argument)
220 else:
221 if not isinstance(c, Indexed):
--> 222 error("Invalid coefficient type for %s" % ufl_err_str(c))
223 f, i = c.ufl_operands
224 if not isinstance(f, Coefficient):
File /usr/local/lib/python3.10/dist-packages/ufl/log.py:158, in Logger.error(self, *message)
156 "Write error message and raise an exception."
157 self._log.error(*message)
--> 158 raise self._exception_type(self._format_raw(*message))
UFLException: Invalid coefficient type for
```
Does anyone know what the problem with `phi` thus defined is?
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [August 4, 2022, 9:38am UTC](https://fenicsproject.discourse.group/t/derive-weak-form-from-variational-principle/8894/6 "2022-08-04T09:38:34Z")
> [@pamgur](#):
>
> `m = ufl.TrialFunction(V)`
`m` should probably be a `dolfinx.fem.Function`
---
### Author:  [@pamgur](https://fenicsproject.discourse.group/u/pamgur)
#### Post date: [August 4, 2022, 4:12pm UTC](https://fenicsproject.discourse.group/t/derive-weak-form-from-variational-principle/8894/7 "2022-08-04T16:12:26Z")
Thanks @dokken for your help! I have more problems now, but I think it’s better to start a new topic.
| | |