Comments (8)
@pierricmora I don’t think we can support the behavior that you are suggesting without adding some very expensive checks to the gmshio readers. One would have to create the intersection of Physical Surfaces tagged by Gmsh.
You can do this by modifiying: https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/io/gmshio.py#L278 to be the unique set of rows in cells.
However, i am not certain the behavior you are suggesting would work for MeshTags of cells, https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/io/gmshio.py#L297
as there is an implicit assumption that a MeshTags object should have a unique set of indices, which is not true in your overlapping case (
dolfinx/cpp/dolfinx/mesh/MeshTags.h
Line 35 in 8e4122b
from dolfinx.
I think that the problem comes from the fact that the cells are read based on the physical groups, rather than by directly asking gmsh to return them. I mean: in extract_topology_and_markers
, the for
loop first gets the dimensions and tags in model.getPhysicalGroups()
, and then only look for cells that share the tag with (entity_types, entity_tags, entity_topologies) = model.mesh.getElements(dim, tag=entity)
. So if a cell is tagged twice, it is returned twice. But also, if a cell is not tagged, then it is never returned.
Could the for loop be reversed, by first reading all cells regardless of tags with model.mesh.getElements(dim, tag=-1)
, and then looking for markers?
from dolfinx.
I think that the problem comes from the fact that the cells are read based on the physical groups, rather than by directly asking gmsh to return them. I mean: in
extract_topology_and_markers
, thefor
loop first gets the dimensions and tags inmodel.getPhysicalGroups()
, and then only look for cells that share the tag with(entity_types, entity_tags, entity_topologies) = model.mesh.getElements(dim, tag=entity)
. So if a cell is tagged twice, it is returned twice. But also, if a cell is not tagged, then it is never returned.Could the for loop be reversed, by first reading all cells regardless of tags with
model.mesh.getElements(dim, tag=-1)
, and then looking for markers?
You can write your own wrapper for this if you would like to. As it is using the GMSH Python API you can create your own custom reader.
Adding all cells regardless of tags is not necessarily a good idea, as GMSH usually only generates cells/facets/edges for Physical entities, to avoid duplicated cells (as in your case).
From section 1.2.3 of https://gmsh.info/doc/texinfo/gmsh.html
1.2.3 Elementary entities vs. physical groups
It is usually convenient to combine elementary geometrical entities into more meaningful groups, e.g. to define some mathematical (“domain”, “boundary with Neumann condition”), functional (“left wing”, “fuselage”) or material (“steel”, “carbon”) properties. Such grouping is done in Gmsh’s geometry module (see Geometry module) through the definition of “physical groups”.
By default in the native Gmsh MSH mesh file format (see Gmsh file formats), as well as in most other mesh formats, if physical groups are defined, the output mesh only contains those elements that belong to at least one physical group. (Different mesh file formats treat physical groups in slightly different ways, depending on their capability to define groups.) To save all mesh elements whether or not physical groups are defined, use the Mesh.SaveAll option (see Mesh options) or specify -save_all on the command line. In some formats (e.g. MSH2), setting Mesh.SaveAll will however discard all physical group definitions.
As I previously stated, our MeshTags
objects will not be able to read in duplicate elements in a sensible way, so rather having a slightly more verbose GMSH tagging routine, as described in
phs1_in_srce = [e for e in phs1 if e in srce] phs1_out_srce = [e for e in phs1 if not(e in srce)] phs2_in_srce = [e for e in phs2 if e in srce] phs2_out_srce = [e for e in phs2 if not(e in srce)] model.add_physical_group(2, phs1_in_srce, 10) model.add_physical_group(2, phs1_out_srce, 11) model.add_physical_group(2, phs2_in_srce, 20) model.add_physical_group(2, phs2_out_srce, 21) model.setPhysicalName(2, 10, "Phase 1 - inside source") model.setPhysicalName(2, 11, "Phase 1 - outside source") model.setPhysicalName(2, 20, "Phase 2 - inside source") model.setPhysicalName(2, 21, "Phase 2 - outside source")
would be preferable from my point of view.
from dolfinx.
Ok, yes it makes sense from Gmsh doc. In my case
e = model.mesh.getElements(dim=2, tag=-1)
print(e[1][0].shape) # (154,): the correct number of cells
does produce the same output regardless the way I have defined - or not - physical groups, so I thought it could be the stable way to go. But indeed, given what's stated in the doc, maybe not, so I think I'll follow your advice. This is a pity, because I find it so much more convenient to work with overlapping or incomplete tags.
from dolfinx.
@jorgensd Then, what about adding these simple checks and error messages to avoid crashes and give the user a hint on what is expected?
The following seems to work on my side. In model_to_mesh
, between lines 258 and 260, add:
# Sort elements by ascending dimension
perm_sort = np.argsort(cell_dimensions)
# Check that all cells are tagged once
_d = model.getDimension()
assert _d in topologies.keys(), 'All cells are expected to be tagged once; none found'
_elementTypes, _elementTags, _nodeTags = model.mesh.getElements(dim=_d, tag=-1)
# assert len(_elementTypes) == 1 # assert only one type of elements; already checked in extract_topology_and_markers
nbcells = len(_elementTags[0])
nbcells_tagged = len(topologies[_d]["topology"])
nbcells_tagged_once = len(np.unique(topologies[_d]['topology'], axis=0))
assert nbcells == nbcells_tagged, f'All cells are expected to be tagged once; found: {nbcells_tagged}, expected: {nbcells}'
assert nbcells_tagged == nbcells_tagged_once, 'All cells are expected to be tagged once; found duplicates'
# Broadcast cell type data and geometric dimension
from dolfinx.
The only thing that worries me here is that the call to nbcells_tagged_once
can be quite expensive, as illustrated in the following mwe:
import numpy as np
from time import perf_counter
N = 1_000_000
M = 8
a = np.arange(N*M).reshape(N, M).astype(np.int64)
b = np.vstack([a, a, a]).astype(np.int64)
start = perf_counter()
unique_b = np.unique(b, axis=0)
end = perf_counter()
print(end-start)
print(f"{unique_b.shape=}, {a.shape=}, {b.shape=}")
yielding
1.6669321240005956
unique_b.shape=(1000000, 8), a.shape=(1000000, 8), b.shape=(3000000, 8)
It would be interesting to see how this modification changes the runtime of "standard" problems. Feel free to make a pull request, and test some standard meshes.
I guess worst case scenario, one has to run Python with -O
to turn of the assert.
from dolfinx.
By playing with N on your snippet I get a linear behaviour for the time spent in np.unique.
Then, I extended my example by adding a mesh_size_max
input parameter for model generation, and a mute_assert
bool in my modified model_to_mesh
to turn on/off the assertions, including np.unique. I ran a few tests representing from 41k to 12M cells (mesh_size_max = 0.01 to 0.01 / 16
) and I obtained a relatively constant +15% increase, unless for the smallest mesh which resulted in a +30% increase (the trend is a slow decrease: largest mesh got a +10%). Is it acceptable? If yes I'll make the PR.
Below is the code I run.
mesh_size_max = 0.01
# Create model
start = perf_counter()
model = gmsh_rect_ellipse(mode=3, mesh_size_max=mesh_size_max, show=False)
end = perf_counter()
gen_model_time = end - start
# model_to_mesh, original
start = perf_counter()
mesh, ct, ft = model_to_mesh(model, MPI.COMM_WORLD, rank=0, gdim=2, mute_assert=True)
end = perf_counter()
model_to_mesh_woutAssert_time = end - start
# model_to_mesh with checks
start = perf_counter()
mesh, ct, ft = model_to_mesh(model, MPI.COMM_WORLD, rank=0, gdim=2, mute_assert=False)
end = perf_counter()
model_to_mesh_withAssert_time = end - start
# say hello
print(f'mesh_size_max: {mesh_size_max},\n' + \
f'number of cells: {len(model.mesh.getElements(dim=2, tag=-1)[1][0])},\n' + \
f'model generation time (Gmsh): {gen_model_time:.3f}s,\n' + \
f'model_to_mesh time, original version: {model_to_mesh_woutAssert_time:.3f}s,\n' + \
f'model_to_mesh time, with checks: {model_to_mesh_withAssert_time:.3f}s,\n' + \
f'overcost due to checks: {(model_to_mesh_withAssert_time - model_to_mesh_woutAssert_time) / model_to_mesh_woutAssert_time * 100:.1f}%')
Outputs:
mesh_size_max: 0.01,
number of cells: 47142,
model generation time (Gmsh): 0.843s,
model_to_mesh time, original version: 0.175s,
model_to_mesh time, with checks: 0.223s,
overcost due to checks: 26.9%
mesh_size_max: 0.005,
number of cells: 186674,
model generation time (Gmsh): 3.652s,
model_to_mesh time, original version: 0.913s,
model_to_mesh time, with checks: 1.035s,
overcost due to checks: 13.4%
mesh_size_max: 0.0025,
number of cells: 743524,
model generation time (Gmsh): 15.690s,
model_to_mesh time, original version: 4.133s,
model_to_mesh time, with checks: 4.747s,
overcost due to checks: 14.9%
mesh_size_max: 0.00125,
number of cells: 2963030,
model generation time (Gmsh): 69.363s,
model_to_mesh time, original version: 21.390s,
model_to_mesh time, with checks: 24.306s,
overcost due to checks: 13.6%
mesh_size_max: 0.000625,
number of cells: 11841660,
model generation time (Gmsh): 329.698s,
model_to_mesh time, original version: 122.245s,
model_to_mesh time, with checks: 134.240s,
overcost due to checks: 9.8%
from dolfinx.
I would be fine with this. 10% overhead isn’t to bad, as one should preferrably only use this wrapper once and not on hpc systems, as Gmsh only writes serial files and Gmsh.model doesn’t do mpi.
from dolfinx.
Related Issues (20)
- Performance regression in collapsing blocked subspace HOT 11
- [BUG]: Visualization of one dimensional data with ParaView HOT 5
- [BUG]: Implicit use of previously created function space HOT 1
- [BUG]: Can't (easily) assemble a block vector without providing a left-hand side in Python HOT 2
- Parent cells not correct in refine_plaza HOT 2
- Access NewtonSolver residual from Python
- Enable macOS CI on PRs (PETSc-less)
- [BUG]: XDMFFile with ASCII encoding is broken in parallel
- Clarity requested: Confusing namespaces `functionspace` is a module in `ufl` but a function in `fem` HOT 6
- [BUG]: dolfinx.io.gmshio.read_from_msh went wrong while using MPI HOT 2
- [BUG]: Unable to run docker container inside a venv HOT 1
- Handling of assertions in python module HOT 1
- Time stamps for XDMF mesh
- Make redhat pipeline compatbile with cffi version 1.17.0
- Lagrange Multiplier for the NS equation with additional constrain
- [BUG]: Calling `create_functionspace` with incompatible element causes breakdown
- oneapi CI failing
- [BUG]: Reading XMDF with prism elements
- Newton solver for blocked problems HOT 2
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from dolfinx.