From 439355ce1c5932bbe2e77ae42e6c00f0c4cdae75 Mon Sep 17 00:00:00 2001
From: Martin Genet <martin.genet@polytechnique.edu>
Date: Tue, 13 Feb 2018 21:20:43 +0100
Subject: [PATCH] getting problem dimension from mesh instead of images
---
fedic.py | 27 +++++++++++++++++----------
1 file changed, 17 insertions(+), 10 deletions(-)
diff --git a/fedic.py b/fedic.py
index 95856c7..15c2238 100644
--- a/fedic.py
+++ b/fedic.py
@@ -112,6 +112,10 @@ def fedic(
assert (os.path.exists(mesh_filename)),\
"No mesh in "+mesh_filename+". Aborting."
mesh = dolfin.Mesh(mesh_filename)
+ mesh_dimension = mesh.ufl_domain().geometric_dimension()
+ assert (mesh_dimension in (2,3)),\
+ "mesh_dimension ("+str(mesh_dimension)+") must be 2 or 3. Aborting."
+ mypy.print_var("mesh_dimension",mesh_dimension,tab+1)
mypy.print_var("mesh_n_cells",len(mesh.cells()),tab+1)
dV = dolfin.Measure(
"dx",
@@ -124,6 +128,9 @@ def fedic(
domain=mesh)
mesh_V0 = dolfin.assemble(dolfin.Constant(1) * dV)
mypy.print_sci("mesh_V0",mesh_V0,tab+1)
+ mesh_h0 = mesh_V0**(1./mesh_dimension)
+ mypy.print_sci("mesh_h0",mesh_h0,tab+1)
+ mesh_h0 = dolfin.Constant(mesh_h0)
if (print_refined_mesh):
mesh_for_plot = dolfin.refine(mesh)
@@ -157,9 +164,9 @@ def fedic(
images_dimension = myvtk.getImageDimensionality(
image=ref_image,
verbose=0)
+ assert (images_dimension == mesh_dimension),\
+ "images_dimension ("+str(images_dimension)+") ≠mesh_dimension ("+str(mesh_dimension)+"). Aborting."
mypy.print_var("images_dimension",images_dimension,tab+1)
- assert (images_dimension in (2,3)),\
- "images_dimension must be 2 or 3. Aborting."
fe = dolfin.FiniteElement(
family="Quadrature",
cell=mesh.ufl_cell(),
@@ -176,7 +183,7 @@ def fedic(
degree=images_quadrature,
quad_scheme="default")
te._quad_scheme = "default" # should not be needed
- for k in xrange(images_dimension**2): # should not be needed
+ for k in xrange(mesh_dimension**2): # should not be needed
te.sub_elements()[k]._quad_scheme = "default" # should not be needed
if (images_expressions_type == "cpp"):
Iref = dolfin.Expression(
@@ -278,7 +285,7 @@ def fedic(
os.remove(filename)
mypy.print_str("Initializing volume file…",tab)
- I = dolfin.Identity(images_dimension)
+ I = dolfin.Identity(mesh_dimension)
F = I + dolfin.grad(U)
J = dolfin.det(F)
@@ -294,12 +301,12 @@ def fedic(
mypy.print_str("Defining regularization energy…",tab)
E = dolfin.Constant(1.0)
nu = dolfin.Constant(regul_poisson)
- kappa = E/3/(1-2*nu) # = E/3 if nu = 0
- lmbda = E*nu/(1+nu)/(1-(images_dimension-1)*nu) # = 0 if nu = 0
- mu = E/2/(1+nu) # = E/2 if nu = 0
- C1 = mu/2 # = E/4 if nu = 0
- C2 = mu/2 # = E/4 if nu = 0
- D1 = kappa/2 # = E/6 if nu = 0
+ kappa = E/3/(1-2*nu) # = E/3 if nu = 0
+ lmbda = E*nu/(1+nu)/(1-(mesh_dimension-1)*nu) # = 0 if nu = 0
+ mu = E/2/(1+nu) # = E/2 if nu = 0
+ C1 = mu/2 # = E/4 if nu = 0
+ C2 = mu/2 # = E/4 if nu = 0
+ D1 = kappa/2 # = E/6 if nu = 0
if (regul_model == "linear"): # <- super bad
e = dolfin.sym(dolfin.grad(U))
--
GitLab