create mesh:
import gmsh
from mpi4py import MPI
gmsh.initialize()
proc = MPI.COMM_WORLD.rank
concrete_marker = 11
rebar_marker = 12
if proc == 0:
gmsh.model.add("concreteslab")
concrete = gmsh.model.occ.add_box(-12, -3, -72, 24, 6, 144)
bary = 2
gmsh.option.setNumber('Mesh.MeshSizeFromCurvature', 5)
cylinder1 = gmsh.model.occ.addCylinder(-10, -bary, 72, 0, 0, -144, .5)
cylinder2 = gmsh.model.occ.addCylinder(0, -bary, 72, 0, 0, -144, .5)
cylinder3 = gmsh.model.occ.addCylinder(10, -bary, 72, 0, 0, -144, .5)
gmsh.model.occ.synchronize()
volTags = gmsh.model.getEntities(3)
gmsh.model.occ.fragment(volTags, volTags)
gmsh.model.addPhysicalGroup(3, [concrete], concrete_marker)
gmsh.model.addPhysicalGroup(3, [cylinder1, cylinder2, cylinder3], rebar_marker)
gmsh.model.mesh.generate(3)
gmsh.write("mesh3D.msh")
gmsh.finalize()
tutorial_linearelasticity:
+*In[ ]:*+
[source, ipython3]
----
# Scaled variable
L = 1
W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma
import numpy as np
import ufl
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from dolfinx import mesh, fem, plot, io
from dolfinx.io import gmshio
domain, cell_markers, facet_markers = gmshio.read_from_msh("mesh3D.msh", MPI.COMM_WORLD, gdim=3)
#domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, W, W])],
# [20,6,6], cell_type=mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("CG", 1))
def clamped_boundary(x):
return np.isclose(x[0], -72)
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = np.array([0,0,0], dtype=ScalarType)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)
T = fem.Constant(domain, ScalarType((0, 0, 0)))
ds = ufl.Measure("ds", domain=domain)
def epsilon(u):
return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*epsilon(u)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, ScalarType((0, 0, -rho*g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
import pyvista
pyvista.start_xvfb()
# Create plotter and pyvista grid
p = pyvista.Plotter()
topology, cell_types, geometry = plot.create_vtk_mesh(V)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
# Attach vector values to grid and warp grid by vector
grid["u"] = uh.x.array.reshape((geometry.shape[0], 3))
actor_0 = p.add_mesh(grid, style="wireframe", color="k")
warped = grid.warp_by_vector("u", factor=1.5)
actor_1 = p.add_mesh(warped, show_edges=True)
p.show_axes()
if not pyvista.OFF_SCREEN:
p.show()
else:
figure_as_array = p.screenshot("deflection.png")
----
+*In[ ]:*+
[source, ipython3]
----
----
+*In[ ]:*+
[source, ipython3]
----
# Scaled variable
L = 1
W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma
import numpy as np
import ufl
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from dolfinx import mesh, fem, plot, io
from dolfinx.io import gmshio
domain, cell_markers, facet_markers = gmshio.read_from_msh("mesh3D.msh", MPI.COMM_WORLD, gdim=3)
#domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, W, W])],
# [20,6,6], cell_type=mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("CG", 1))
def clamped_boundary(x):
return np.isclose(x[0], -72)
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = np.array([0,0,0], dtype=ScalarType)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)
T = fem.Constant(domain, ScalarType((0, 0, 0)))
ds = ufl.Measure("ds", domain=domain)
def epsilon(u):
return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*epsilon(u)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, ScalarType((0, 0, -rho*g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
import pyvista
pyvista.start_xvfb()
# Create plotter and pyvista grid
p = pyvista.Plotter()
topology, cell_types, geometry = plot.create_vtk_mesh(V)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
# Attach vector values to grid and warp grid by vector
grid["u"] = uh.x.array.reshape((geometry.shape[0], 3))
actor_0 = p.add_mesh(grid, style="wireframe", color="k")
warped = grid.warp_by_vector("u", factor=1.5)
actor_1 = p.add_mesh(warped, show_edges=True)
p.show_axes()
if not pyvista.OFF_SCREEN:
p.show()
else:
figure_as_array = p.screenshot("deflection.png")
----
ValueError:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
File /usr/local/lib/python3.10/dist-packages/IPython/core/formatters.py:972, in MimeBundleFormatter.__call__(self, obj, include, exclude)
969 method = get_real_method(obj, self.print_method)
971 if method is not None:
--> 972 return method(include=include, exclude=exclude)
973 return None
974 else:
File /usr/lib/python3/dist-packages/panel/viewable.py:692, in Viewable._repr_mimebundle_(self, include, exclude)
690 if config.embed:
691 return render_model(model)
--> 692 return render_mimebundle(model, doc, comm, manager, location)
File /usr/lib/python3/dist-packages/panel/io/notebook.py:197, in render_mimebundle(model, doc, comm, manager, location)
195 loc = location._get_model(doc, model, model, comm)
196 doc.add_root(loc)
--> 197 return render_model(model, comm)
File /usr/lib/python3/dist-packages/panel/io/notebook.py:160, in render_model(model, comm)
156 from ..config import panel_extension as pnext
158 target = model.ref['id']
--> 160 (docs_json, [render_item]) = standalone_docs_json_and_render_items([model], True)
161 div = div_for_render_item(render_item)
162 render_json = render_item.to_json()
File /usr/lib/python3/dist-packages/bokeh/embed/util.py:337, in standalone_docs_json_and_render_items(models, suppress_callback_warning)
335 docs_json: Dict[ID, DocJson] = {}
336 for doc, (docid, _) in docs.items():
--> 337 docs_json[docid] = doc.to_json()
339 render_items: List[RenderItem] = []
340 for _, (docid, roots) in docs.items():
File /usr/lib/python3/dist-packages/bokeh/document/document.py:758, in Document.to_json(self)
749 ''' Convert this document to a JSON object.
750
751 Return:
752 JSON-data
753
754 '''
756 # this is a total hack to go via a string, needed because
757 # our BokehJSONEncoder goes straight to a string.
--> 758 doc_json = self.to_json_string()
759 return loads(doc_json)
File /usr/lib/python3/dist-packages/bokeh/document/document.py:790, in Document.to_json_string(self, indent)
778 root_ids = [ r.id for r in self._roots ]
780 json = DocJson(
781 title=self.title,
782 defs=serializer.definitions,
(...)
787 version=__version__,
788 )
--> 790 return serialize_json(json, indent=indent)
File /usr/lib/python3/dist-packages/bokeh/core/json_encoder.py:171, in serialize_json(obj, pretty, indent, **kwargs)
168 if pretty and indent is None:
169 indent = 2
--> 171 return json.dumps(obj, cls=BokehJSONEncoder, allow_nan=False, indent=indent, separators=separators, sort_keys=True, **kwargs)
File /usr/lib/python3.10/json/__init__.py:238, in dumps(obj, skipkeys, ensure_ascii, check_circular, allow_nan, cls, indent, separators, default, sort_keys, **kw)
232 if cls is None:
233 cls = JSONEncoder
234 return cls(
235 skipkeys=skipkeys, ensure_ascii=ensure_ascii,
236 check_circular=check_circular, allow_nan=allow_nan, indent=indent,
237 separators=separators, default=default, sort_keys=sort_keys,
--> 238 **kw).encode(obj)
File /usr/lib/python3.10/json/encoder.py:199, in JSONEncoder.encode(self, o)
195 return encode_basestring(o)
196 # This doesn't pass the iterator directly to ''.join() because the
197 # exceptions aren't as detailed. The list call should be roughly
198 # equivalent to the PySequence_Fast that ''.join() would do.
--> 199 chunks = self.iterencode(o, _one_shot=True)
200 if not isinstance(chunks, (list, tuple)):
201 chunks = list(chunks)
File /usr/lib/python3.10/json/encoder.py:257, in JSONEncoder.iterencode(self, o, _one_shot)
252 else:
253 _iterencode = _make_iterencode(
254 markers, self.default, _encoder, self.indent, floatstr,
255 self.key_separator, self.item_separator, self.sort_keys,
256 self.skipkeys, _one_shot)
--> 257 return _iterencode(o, 0)
ValueError: Out of range float values are not JSON compliant
There is a value error:
actor_1 = p.add_mesh(warped, show_edges=True)
What I am trying to do is get the newly create and imported mesh3D working the same way as does the fenicsx tutorial example. Is there a way to resolve the ValueError that is coming out and get a warped mesh the same way for mesh3D?