From a826939f6638bdd6cbc7a84b22dfba0eed119c8e Mon Sep 17 00:00:00 2001
From: karamoko <karamoko.samassa@inria.fr>
Date: Fri, 14 Oct 2022 11:48:19 +0200
Subject: [PATCH] adding variable dirichlet with user expression
---
src/bvpy/boundary_conditions/dirichlet.py | 30 ++++++++++++---------
tutorials/bvpy_tutorial_1_hello_world.ipynb | 23 ++++++++++++----
2 files changed, 36 insertions(+), 17 deletions(-)
diff --git a/src/bvpy/boundary_conditions/dirichlet.py b/src/bvpy/boundary_conditions/dirichlet.py
index 0ab15c1..22fdca0 100644
--- a/src/bvpy/boundary_conditions/dirichlet.py
+++ b/src/bvpy/boundary_conditions/dirichlet.py
@@ -24,15 +24,15 @@ import numpy as np
from bvpy import logger
from dolfinx import fem, mesh
from petsc4py.PETSc import ScalarType
-# from sympy import Symbol
-# from sympy.parsing.sympy_parser import parse_expr
+from sympy import Symbol, symbols, lambdify
+from sympy.parsing.sympy_parser import parse_expr
from .boundary import Boundary
from bvpy.domains.geometry import boundary_normal
-# _XYZ = {'x': Symbol('x[0]'),
-# 'y': Symbol('x[1]'),
-# 'z': Symbol('x[2]')}
+_XYZ = {'x': Symbol('x[0]'),
+ 'y': Symbol('x[1]'),
+ 'z': Symbol('x[2]')}
class ConstantDirichlet(object):
@@ -187,7 +187,7 @@ class VariableDirichlet(object):
The boundary condition in a fenics-readable format.
"""
- # parse = parse_expr(self._val, local_dict=_XYZ, evaluate=False)
+ parse = parse_expr(self._val, local_dict=_XYZ, evaluate=False)
if self._degree is None:
# expr = fe.Expression(str(parse),
# element=functionSpace.ufl_element())
@@ -201,13 +201,19 @@ class VariableDirichlet(object):
# self._boundary(),
# method=self._method)
- # Need this one here
- # uD = fem.Function(functionSpace)
- # uD.interpolate(lambda x: 1 + x[0] - x[1])
+ _x = symbols("x")
+ _lam_f = lambdify(_x, parse)
+
+ u_dir = fem.Function(functionSpace)
+ u_dir.interpolate(_lam_f)
- boundary_dofs = fem.locate_dofs_geometrical(functionSpace, self._boundary())
- dir = fem.dirichletbc(ScalarType(0), boundary_dofs, functionSpace)
- print(boundary_dofs)
+ _domain = functionSpace.mesh
+ tdim = _domain.topology.dim
+ _domain.topology.create_connectivity(tdim-1, tdim)
+ boundary_factes = mesh.exterior_facet_indices(_domain.topology)
+ boundary_dofs = fem.locate_dofs_topological(functionSpace, tdim- 1, boundary_factes)
+
+ dir = fem.dirichletbc(u_dir, boundary_dofs)
return dir
else:
# dir = fe.DirichletBC(functionSpace.sub(self._subspace),
diff --git a/tutorials/bvpy_tutorial_1_hello_world.ipynb b/tutorials/bvpy_tutorial_1_hello_world.ipynb
index 45e1c4a..f0268e5 100644
--- a/tutorials/bvpy_tutorial_1_hello_world.ipynb
+++ b/tutorials/bvpy_tutorial_1_hello_world.ipynb
@@ -119,8 +119,8 @@
"outputs": [],
"source": [
"from bvpy.boundary_conditions import dirichlet\n",
- "# diri = dirichlet('1 + x*x + 2*y*y', boundary='all')\n",
- "diri = dirichlet(-5, boundary='all')"
+ "diri = dirichlet('1 + x*x + 2*y*y', boundary='all')\n",
+ "# diri = dirichlet(-5, boundary='all')"
]
},
{
@@ -141,7 +141,20 @@
"cell_type": "code",
"execution_count": 4,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[ 0 2 3 6 10 11 15 18 23 27 32 33 38 45 46 52 54 63\n",
+ " 64 81 82 91 92 104 117 118 120 122 124 126 133 150 181 182 208 223\n",
+ " 236 249 259 279 286 307 315 329 341 359 383 388 410 415 437 444 445 446\n",
+ " 456 457 458 468 478 480 486 493 495 509 512 515 516 523 528 532 533 534\n",
+ " 537 542 544 546 549 550 551 553]\n",
+ "<dolfinx.fem.bcs.DirichletBC object at 0x190b8d1c0>\n"
+ ]
+ }
+ ],
"source": [
"from bvpy import BVP\n",
"prblm_1 = BVP(domain= rectangle, vform=poisson, bc=diri)"
@@ -197,7 +210,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
- "model_id": "b954cfda580d4d3e81eb3155260ff58a",
+ "model_id": "5d41246f0b9e475a853df70faa45a850",
"version_major": 2,
"version_minor": 0
},
@@ -251,7 +264,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
- "Maximal discrepency between expectation and computed solution: 9.00e+00\n"
+ "Maximal discrepency between expectation and computed solution: 5.16e-01\n"
]
}
],
--
GitLab