diff --git a/src/bvpy/bvp.py b/src/bvpy/bvp.py
index 2f3299b77c79255052d49170916ac0986d06f7aa..246b53bef2784daa72501e482f34377a7899cfcd 100644
--- a/src/bvpy/bvp.py
+++ b/src/bvpy/bvp.py
@@ -456,14 +456,22 @@ class BVP(object):
# self.solver.set_problem(self.solution, lhs,
# rhs, self.dirichletCondition,
# self._pointsources)
-
- self.solver = fem.petsc.LinearProblem(lhs,
- rhs,
- bcs=self.dirichletCondition,
- petsc_options={
- "ksp_type": "preonly",
- "pc_type": "lu"})
- self.solution = self.solver.solve()
+ if isinstance(rhs, ufl.constantvalue.Zero):
+ self.solver = NonLinearSolver()
+ self._set_solver_parameter(**kwargs)
+ jacobian = ufl.derivative(lhs, self.solution, self.trialFunction)
+ self.solver.set_problem(self.solution, lhs,
+ jacobian, self.dirichletCondition,
+ self._pointsources)
+ self.solver.solve()
+ else:
+ self.solver = fem.petsc.LinearProblem(lhs,
+ rhs,
+ bcs=self.dirichletCondition,
+ petsc_options={
+ "ksp_type": "preonly",
+ "pc_type": "lu"})
+ self.solution = self.solver.solve()
def info(self):
"""Prints the descrition of the problem.
diff --git a/src/bvpy/utils/post_processing.py b/src/bvpy/utils/post_processing.py
index 6a0babe9393227785a0c6fa4c2b6e84ff4395a68..95cdcd82632263970cfa13b37724db6b8fe70643 100644
--- a/src/bvpy/utils/post_processing.py
+++ b/src/bvpy/utils/post_processing.py
@@ -25,6 +25,7 @@ import numpy as np
import operator
from bvpy.utils.pre_processing import create_expression
+import basix.ufl_wrapper as bu
import ufl
from dolfinx import fem
@@ -75,14 +76,15 @@ class SolutionExplorer(object):
if self._dof_per_nodes == 0:
self._dof_per_nodes = int(1) # for scalar Function
- self.nb_dofs = len(self._functionSpace.dofmap().dofs())
+ # self.nb_dofs = len(self._functionSpace.dofmap().dofs())
+ self.nb_dofs = len(self._functionSpace.dofmap.list)
self.nb_nodes = int(self.nb_dofs/self._dof_per_nodes)
- if isinstance(self._functionSpace.ufl_element(), ufl.FiniteElement):
+ if isinstance(self._functionSpace.ufl_element(), bu.BasixElement):
self._type = 'scalar'
- elif isinstance(self._functionSpace.ufl_element(), ufl.VectorElement):
+ elif isinstance(self._functionSpace.ufl_element(), bu.VectorElement):
self._type = 'vector'
- elif isinstance(self._functionSpace.ufl_element(), ufl.TensorElement):
+ elif isinstance(self._functionSpace.ufl_element(), bu.TensorElement):
self._type = 'tensor'
def __call__(self, *x):
@@ -97,7 +99,8 @@ class SolutionExplorer(object):
Array containing the vertex position
"""
- return self._mesh.coordinates()
+ # return self._mesh.coordinates()
+ return self._mesh.geometry.x
def get_dofs_position(self):
"""Generates a Fenics object containing the dofs positions.
@@ -119,7 +122,8 @@ class SolutionExplorer(object):
The array of the dofs' values
"""
- return self._fenics.vector().get_local()
+ # return self._fenics.vector().get_local()
+ return self._fenics.vector.array
def to_fenics(self):
"""Returns a fenics.Function from the considered function.
@@ -141,23 +145,31 @@ class SolutionExplorer(object):
"""
if self._type == 'scalar':
- return self._fenics.compute_vertex_values()
+ # return self._fenics.compute_vertex_values()
+ return self._fenics.x.array
else:
size_n = self.shape[0]
if self._type == 'vector':
- res = np.zeros((self._mesh.num_entities(0), self.shape[0]))
+ # res = np.zeros((self._mesh.num_entities(0), self.shape[0]))
+ res = np.zeros((self._mesh.geometry.x.shape[0], self.shape[0]))
for i in range(size_n):
- res[:, i] = self._fenics.sub(i, deepcopy=True).compute_vertex_values()
+ # res[:, i] = self._fenics.sub(i, deepcopy=True).compute_vertex_values()
+ res[:, i] = self._fenics.sub(i).collapse().x.array
if self._type == 'tensor':
- res = np.zeros((self._mesh.num_entities(0),
+ # res = np.zeros((self._mesh.num_entities(0),
+ # self.shape[0], self.shape[1]))
+ res = np.zeros((self._fenics.x.array.shape[0],
self.shape[0], self.shape[1]))
count = 0
for i in range(size_n):
for j in range(size_n):
- res[:, i, j] = self._fenics.sub(count+j, deepcopy=True).compute_vertex_values()
+ # res[:, i, j] = self._fenics.sub(count+j, deepcopy=True).compute_vertex_values()
+ #TODO: BE CAREFUL HERE j is no more used
+ res[:, i, j] = self._fenics.x.array
count += size_n
return res
+
def geometric_filter(self, axe='x', threshold=0, comparator=">",
return_position=False, sort=None):
diff --git a/src/bvpy/vforms/elasticity.py b/src/bvpy/vforms/elasticity.py
index e02110e68551ef82dcea55aaa7427ca63b459b1e..66d92199b6f3e84e488ba9fe25d027c7e2c7aea2 100644
--- a/src/bvpy/vforms/elasticity.py
+++ b/src/bvpy/vforms/elasticity.py
@@ -329,10 +329,11 @@ class HyperElasticForm(ElasticForm):
def construct_form(self, u, v, sol):
E = self._strain(sol)
+ mesh = v.ufl_function_space().mesh
# -- Constitutive models
if self.model == 'st_venant':
- W = ufl.Constant(0.5) * self._parameters['lmbda']\
+ W = fem.Constant(mesh, 0.5) * self._parameters['lmbda']\
* ufl.tr(E) * ufl.tr(E)
W += self._parameters['mu'] * ufl.tr(E * E)
@@ -349,7 +350,12 @@ class HyperElasticForm(ElasticForm):
ufl.ln(J) + (self._parameters['lmbda'] / 2) * (ufl.ln(J))**2
self.lhs = self._parameters['thickness'] * W * self.dx
- self.lhs -= ufl.dot(self._parameters['source'], sol) * self.dx
+
+ source = self._parameters['source']
+ if isinstance(source, (int, float, list, tuple, np.ndarray)):
+ source = fem.Constant(mesh, source)
+
+ self.lhs -= ufl.dot(source, sol) * self.dx
self.lhs = ufl.derivative(self.lhs, sol, v)
diff --git a/tutorials/bvpy_tutorial_7_hyper_elasticity.ipynb b/tutorials/bvpy_tutorial_7_hyper_elasticity.ipynb
index 1d8ed92fac62b38c13966ad876dc13e38821f40d..321e2eecfd2ac5c0f1545e5f7b386b07c58c11d2 100644
--- a/tutorials/bvpy_tutorial_7_hyper_elasticity.ipynb
+++ b/tutorials/bvpy_tutorial_7_hyper_elasticity.ipynb
@@ -217,7 +217,7 @@
"source": [
"from bvpy.vforms import LinearElasticForm\n",
"\n",
- "linear_response = LinearElasticForm(young=Young, poisson=Poisson, source=[0., 0.], plane_stress=True)\n",
+ "linear_response = LinearElasticForm(young=Young, poisson=Poisson, source=[0., 0., 0.], plane_stress=True)\n",
"\n",
"linear_response.info()"
]
@@ -237,7 +237,7 @@
"source": [
"from bvpy.vforms import HyperElasticForm\n",
"\n",
- "non_linear_response = HyperElasticForm(young=Young, poisson=Poisson, source=[0., 0.], plane_stress=True)\n",
+ "non_linear_response = HyperElasticForm(young=Young, poisson=Poisson, source=[0., 0., 0.], plane_stress=True)\n",
"\n",
"non_linear_response.info()"
]
@@ -299,10 +299,9 @@
" for name, response in {'ln': linear_response, 'nl': non_linear_response}.items():\n",
" # -- Setting the problem\n",
" stretching = BVP(patch, response, traction)\n",
- " \n",
+ " \n",
" # -- Solving the problems\n",
" stretching.solve()\n",
- " \n",
" # -- Computing the corresponding strain field\n",
" fe_strain = response.strain(stretching.solution)\n",
" np_strain = SolutionExplorer(fe_strain).to_vertex_values()\n",
@@ -463,7 +462,9 @@
},
"outputs": [],
"source": [
- "import fenics as fe\n",
+ "# import fenics as fe\n",
+ "from dolfinx import fem, mesh\n",
+ "\n",
"\n",
"def compute_young(concentrations, Young_iso=1e3, domain=globe, component=0):\n",
" \"\"\"Computes Young's modulus given a concentration field.\n",
@@ -489,23 +490,38 @@
" \n",
" \"\"\"\n",
" \n",
- " vct_func_space = fe.VectorFunctionSpace(domain.mesh, \"P\", 1, dim=2)\n",
+ " # vct_func_space = fe.VectorFunctionSpace(domain.mesh, \"P\", 1, dim=2)\n",
+ "\n",
+ " vct_func_space = fem.VectorFunctionSpace(domain.mesh, (\"P\", 1))\n",
+ " c_func = fem.Function(vct_func_space)\n",
+ "\n",
+ " if isinstance(concentrations, (list, np.ndarray)):\n",
+ " # c_func = fe.interpolate(concentrations, vct_func_space).sub(component)\n",
+ " if(vct_func_space.num_sub_spaces > 0):\n",
+ " for i in range(vct_func_space.num_sub_spaces):\n",
+ " # self.solution.sub(i).interpolate(solution[i])\n",
+ " c_func.sub(i).interpolate(concentrations[i])\n",
+ " else:\n",
+ " c_func.interpolate(concentrations)\n",
" \n",
- " if isinstance(concentrations, fe.UserExpression):\n",
- " c_func = fe.interpolate(concentrations, vct_func_space).sub(component)\n",
- " \n",
- " elif isinstance(concentrations, fe.Function):\n",
- " c_func = fe.project(concentrations, vct_func_space).sub(component)\n",
+ " elif isinstance(concentrations, fem.Function):\n",
+ " # c_func = fe.project(concentrations, vct_func_space).sub(component)\n",
+ " c_func.interpolate(concentrations)\n",
" \n",
" else:\n",
" print(\"WARNING: the type of the input is not supported\")\n",
" \n",
- " c_max = c_func.vector().max()\n",
- " c_min = c_func.vector().min()\n",
- "\n",
- " Young = fe.Constant(1)\n",
- " Young += fe.Constant(100 / (c_max - c_min)) * (fe.Constant(c_max) - c_func)\n",
- " Young *= fe.Constant(Young_iso)\n",
+ " # c_max = c_func.vector().max()\n",
+ " # c_min = c_func.vector().min()\n",
+ " c_max = c_func.vector.max()\n",
+ " c_min = c_func.vector.min()\n",
+ "\n",
+ " # Young = fe.Constant(1)\n",
+ " # Young += fe.Constant(100 / (c_max - c_min)) * (fe.Constant(c_max) - c_func)\n",
+ " # Young *= fe.Constant(Young_iso)\n",
+ " Young = fem.Constant(domain.mesh, 1)\n",
+ " Young += fem.Constant(domain.mesh, 100 / (c_max - c_min)) * (fem.Constant(domain.mesh, c_max) - c_func)\n",
+ " Young *= fem.Constant(domain.mesh, Young_iso)\n",
" \n",
" return Young"
]
@@ -718,19 +734,30 @@
"metadata": {},
"outputs": [],
"source": [
- "from bvpy.utils.pre_processing import create_expression\n",
+ "# from bvpy.utils.pre_processing import create_expression\n",
+ "\n",
+ "# # -- Creating the initial concentration distribution\n",
+ "# ka = 9\n",
+ "# kb = 11\n",
+ "\n",
+ "# np.random.seed(1093)\n",
+ "# globe.discretize()\n",
+ "# random = np.random.uniform(low=-1, high=1, size=globe.mesh.num_entities(2))\n",
+ "\n",
+ "# random_pattern = [[1+0.04*ka**2+0.1*r, 0.2*kb+0.1*r] for r in random]\n",
+ "\n",
+ "# concentrations = create_expression(random_pattern)\n",
+ "\n",
+ "import numpy as np\n",
"\n",
- "# -- Creating the initial concentration distribution\n",
"ka = 9\n",
"kb = 11\n",
"\n",
- "np.random.seed(1093)\n",
"globe.discretize()\n",
- "random = np.random.uniform(low=-1, high=1, size=globe.mesh.num_entities(2))\n",
- "\n",
- "random_pattern = [[1+0.04*ka**2+0.1*r, 0.2*kb+0.1*r] for r in random]\n",
- "\n",
- "concentrations = create_expression(random_pattern)"
+ "dim = 3 # each node is presented by (x, y, z)\n",
+ "_size = globe.mesh.geometry.dofmap.num_nodes*dim\n",
+ "random = np.random.uniform(low=-1, high=1, size=_size)\n",
+ "concentrations = [1+ 0.04*ka**2+0.1*random, 0.2*kb+0.1*random] "
]
},
{