From 6b228690be37ff48df319df23308b020e1f34302 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20G=C3=BCnther?= Date: Fri, 12 Aug 2022 12:15:01 +0200 Subject: [PATCH 1/8] minors --- pygimli/meshtools/polytools.py | 12 ++++++------ pygimli/viewer/showmesh.py | 5 +++-- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/pygimli/meshtools/polytools.py b/pygimli/meshtools/polytools.py index 424a30bb2..20a9b4f67 100644 --- a/pygimli/meshtools/polytools.py +++ b/pygimli/meshtools/polytools.py @@ -1215,10 +1215,10 @@ def createParaMeshPLC3D(sensors, paraDX=0, paraDepth=-1, paraBoundary=None, Margin for parameter domain in relative extend. paraMaxCellSize: double, optional - Maximum cell size for parametric region in m² + Maximum cell size for parametric region in m³ boundaryMaxCellSize: double, optional - Maximum cells size in the boundary region in m² + Maximum cells size in the boundary region in m³ boundary: [float, float] [10., 10.] Boundary width to be appended for domain prolongation in relative @@ -1256,7 +1256,7 @@ def createParaMeshPLC3D(sensors, paraDX=0, paraDepth=-1, paraBoundary=None, paraDepth = np.median(sensors[:, 2]) - paraDepth depth = paraDepth - max(boundary[0]*xSpan, boundary[1]*ySpan)/2 - + def sortP(p): base = pg.core.Line(p[0], p[1]).at(-1e7) def cmp_(p1, p2): @@ -1264,7 +1264,7 @@ def cmp_(p1, p2): return -1 else: return 1 - + p.sort(key=functools.cmp_to_key(cmp_)) bounds = [surface] @@ -1330,8 +1330,8 @@ def cmp_(p1, p2): if paraDX > 0: for s in sensors: pdPLC.createNode(s - [0.0, 0.0, paraDX]) - - pdPLC.addRegionMarker(pg.center(bttmA.positions()) + [0.0, 0.0, 0.1], + + pdPLC.addRegionMarker(pg.center(bttmA.positions()) + [0.0, 0.0, 0.1], marker=1, area=boundaryMaxCellSize) pdPLC.addRegionMarker(pg.center(bttmP.positions()) + [0.0, 0.0, 0.1], marker=2, area=paraMaxCellSize) diff --git a/pygimli/viewer/showmesh.py b/pygimli/viewer/showmesh.py index 0be8b22ec..ed555c669 100644 --- a/pygimli/viewer/showmesh.py +++ b/pygimli/viewer/showmesh.py @@ -328,7 +328,8 @@ def showMesh(mesh, data=None, block=False, colorBar=None, data = pg.solver.parseMapToCellArray(data, mesh) if hasattr(data[0], '__len__') and not \ - isinstance(data, np.ma.core.MaskedArray): + isinstance(data, np.ma.core.MaskedArray) and not \ + isinstance(data[0], str): # [u,v] x N if len(data) == 2: @@ -342,7 +343,7 @@ def showMesh(mesh, data=None, block=False, colorBar=None, elif data.shape[1] == 3: # if sum(data[:, 0]) != sum(data[:, 1]): # drawStreams(ax, mesh, data, **kwargs) - drawStreams(ax, mesh, data[:, 0:2], **kwargs) + drawStreams(ax, mesh, data[:, :2], **kwargs) else: # Try animation frames x N From c8075e8eb8cab71dfdd1b30a19582878abf61cad Mon Sep 17 00:00:00 2001 From: carsten-forty2 Date: Fri, 12 Aug 2022 14:25:13 +0200 Subject: [PATCH 2/8] Add line_width to pv show options, fix example plot_05_refraction_3d.py --- doc/examples/2_seismics/plot_05_refraction_3D.py | 8 +++++--- pygimli/viewer/pv/draw.py | 4 +++- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/doc/examples/2_seismics/plot_05_refraction_3D.py b/doc/examples/2_seismics/plot_05_refraction_3D.py index e81ac3404..549f4b15e 100644 --- a/doc/examples/2_seismics/plot_05_refraction_3D.py +++ b/doc/examples/2_seismics/plot_05_refraction_3D.py @@ -48,7 +48,7 @@ if pyvista: label = pg.utils.unit("vel") - pg.show(mesh, vel, label=label, notebook=True) + pg.show(mesh, vel, label=label, showMesh=True) ################################################################################ # Set-up data container. @@ -81,7 +81,8 @@ # Plot final ray paths. if pyvista: - plotter, _ = pg.show(mesh, hold=True) + plotter, _ = pg.show(mesh, style='wireframe', line_width=0.1, + hold=True) drawSensors(plotter, sensors, diam=0.5, color='yellow') for ray in rays: @@ -90,4 +91,5 @@ stop = tuple(ray[i + 1]) line = pyvista.Line(start, stop) plotter.add_mesh(line, color='green', line_width=2) - # plotter.show() + + plotter.show() diff --git a/pygimli/viewer/pv/draw.py b/pygimli/viewer/pv/draw.py index f53292e2a..9fd340399 100644 --- a/pygimli/viewer/pv/draw.py +++ b/pygimli/viewer/pv/draw.py @@ -42,9 +42,10 @@ def drawMesh(ax, mesh, notebook=False, **kwargs): returnActor = kwargs.pop('returnActor', False) showMesh = kwargs.pop('showMesh', False) grid = kwargs.pop('grid', False) - colorBar = kwargs.pop('colorBar', True) + colorBar = kwargs.pop('colorBar', style != 'wireframe') name = kwargs.pop('name', 'Mesh') bc = kwargs.pop('bc', '#EEEEEE') # background color + lw = kwargs.pop('line_width', 0.1) filt = kwargs.pop('filter', {}) dataName = kwargs.pop('label', list(mesh.cell_data.keys())[0]) @@ -87,6 +88,7 @@ def drawMesh(ax, mesh, notebook=False, **kwargs): color=color, style=style, show_edges=showMesh, + line_width=lw, #edge_color='white', show_scalar_bar=colorBar, opacity=opacity, From c3ecd6be8fcf885b067438f2005fe43cb19376c9 Mon Sep 17 00:00:00 2001 From: carsten-forty2 Date: Fri, 12 Aug 2022 15:29:25 +0200 Subject: [PATCH 3/8] improve Tank3d example a bit --- doc/examples/3_dc_and_ip/plot_modTank3d.py | 23 +++++++++++----------- pygimli/viewer/pv/vistaview.py | 12 +++++------ 2 files changed, 18 insertions(+), 17 deletions(-) diff --git a/doc/examples/3_dc_and_ip/plot_modTank3d.py b/doc/examples/3_dc_and_ip/plot_modTank3d.py index f69d08fe1..462eb7383 100644 --- a/doc/examples/3_dc_and_ip/plot_modTank3d.py +++ b/doc/examples/3_dc_and_ip/plot_modTank3d.py @@ -66,19 +66,21 @@ ############################################################################### # For sufficient numerical accuracy it is generally a good idea to refine the # mesh in the vicinity of the electrodes positions. -# We force the local mesh refinement by an additional node at 1/2 mm +# We force the local mesh refinement by an additional node at 1 mm # distance in -z-direction. for s in plc.positions(pg.find(plc.nodeMarkers() == -99)): - plc.createNode(s - [0.0, 0.0, 1e-3/2]) + plc.createNode(s - [0.0, 0.0, 1e-3]) # Also refine the reference node -plc.createNode([0.5, 0.5, -0.5 - 1e-3/2]) +plc.createNode([0.5, 0.5, -0.5 - 1e-3]) +#pg.show(plc, markers=True, showMesh=True) ############################################################################### # Create the tetrahedron mesh (calling the tetgen mesh generator) mesh = mt.createMesh(plc) +pg.show(mesh, markers=True, showMesh=True) ############################################################################### # First we want to simulate our ERT response for a homogeneous resistivity @@ -114,17 +116,18 @@ cube = mt.createCube(size=[0.3, 0.2, 0.8], pos=[0.7, 0.2], marker=2) plc += cube +pg.show(plc, alpha=0.3) mesh = mt.createMesh(plc) -print(mesh) ############################################################################### -# Until we find a way (TODO) showing 3D structures in a 2D documentation, -# it is advisable to control the geometry in a 3D viewer. We recommend Paraview -# https://www.paraview.org/ and will export the geometry and in the -# vtk file format. +# Also its advisable to control the geometry in a 3D viewer. +# We recommend Paraview https://www.paraview.org/ and +# will export the geometry and in the vtk file format. # plc.exportVTK('plc') # mesh.exportVTK('mesh') +pg.show(mesh, mesh.cellMarkers(), showMesh=True, + filter={'clip':{'origin':(0.7, 0, 0.0)},}) ############################################################################### # Now that we have a mesh with different regions and cell markers, we can define @@ -135,7 +138,6 @@ res = [[1, 10.0], [2, 100.0]] # map markers 1 and 2 to 10 and 100 Ohmm, resp. het = ert.simulate(mesh, res=res, scheme=shm, sr=False, calcOnly=True, verbose=True) -pg.show(mesh, notebook=True) ############################################################################### # The apparent resistivity for a homogeneous model of 1 Ohmm should be 1 Ohmm. @@ -152,8 +154,6 @@ # np.testing.assert_approx_equal(het('k')[0], 0.820615269548) -pg.wait() - ############################################################################### # For such kind of simulations, the homogeneous part should be high accurate # because it is usually needed once after storing the geometric factors. @@ -168,5 +168,6 @@ # checks: # TODO: # * any idea for a Figure here? +# * maybe show data and effect of topography # * Alternatively, we can create a real cavity by changing the marker # in isHole flag for createCube (check) diff --git a/pygimli/viewer/pv/vistaview.py b/pygimli/viewer/pv/vistaview.py index 59af56e79..6a6185548 100644 --- a/pygimli/viewer/pv/vistaview.py +++ b/pygimli/viewer/pv/vistaview.py @@ -117,17 +117,17 @@ def showMesh3DVista(mesh, data=None, **kwargs): if kwargs.get('aa', False): plotter.enable_anti_aliasing() - if notebook: + if notebook is True: # monkey patch show of this plotter instance so we can use multiple backends and only plotter.show() .. whoever this needs. plotter.__show = plotter.show - plotter.show = lambda *args, **kwargs: plotter.__show( - *args, jupyter_backend=backend, **kwargs - ) + plotter.show = lambda *args, **kwargs: plotter.__show(*args, + jupyter_backend=backend, **kwargs) else: plotter.__show = plotter.show - plotter.show = lambda *args, **kwargs: plotter.__show(*args, **kwargs) if pg.viewer.mpl.isInteractive() else False + plotter.show = lambda *args, **kwargs: plotter.__show(*args, + **kwargs) if pg.viewer.mpl.isInteractive() else False - if not hold: + if hold is False: plotter.show() # , None to keep compatability From 2d346c83550456f45bdecd84fb99873ca2b458b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20G=C3=BCnther?= Date: Fri, 12 Aug 2022 18:06:42 +0200 Subject: [PATCH 4/8] change ERT simulate critical to warn for Jenkins --- pygimli/physics/ert/ertManager.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygimli/physics/ert/ertManager.py b/pygimli/physics/ert/ertManager.py index 4d202ae57..726992582 100644 --- a/pygimli/physics/ert/ertManager.py +++ b/pygimli/physics/ert/ertManager.py @@ -200,7 +200,7 @@ def simulate(self, mesh, scheme, res, **kwargs): # >>> rhoa = data.get('rhoa').array() # >>> phia = data.get('phia').array() """ - pg.critical('Obsolete, do not use!. Use ert.simulate instead') + pg.warn('Obsolete, do not use!. Use ert.simulate instead') verbose = kwargs.pop('verbose', self.verbose) calcOnly = kwargs.pop('calcOnly', False) From aac517523a09528cb152f33845c9571ac4720677 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20G=C3=BCnther?= Date: Fri, 12 Aug 2022 18:46:42 +0200 Subject: [PATCH 5/8] clean petro example --- doc/examples/5_misc/plot_petro_joint_inv.py | 38 ++++++++++++--------- pygimli/frameworks/methodManager.py | 22 +++++++----- 2 files changed, 35 insertions(+), 25 deletions(-) diff --git a/doc/examples/5_misc/plot_petro_joint_inv.py b/doc/examples/5_misc/plot_petro_joint_inv.py index df5917f1b..d13396202 100644 --- a/doc/examples/5_misc/plot_petro_joint_inv.py +++ b/doc/examples/5_misc/plot_petro_joint_inv.py @@ -5,8 +5,8 @@ ----------------------------- Joint inversion of different geophysical techniques helps to improve both -resolution and interpretability of the resulting images. Different data sets can -be directly coupled, if there is a link to an underlying target parameter. +resolution and interpretability of the resulting images. Different data sets +can be directly coupled, if there is a link to an underlying target parameter. In this example, ERT and traveltime data are inverted for water saturation. For details see section 3.3 of the pyGIMLi paper (https://cg17.pygimli.org). """ @@ -16,12 +16,15 @@ import pygimli as pg from pygimli import meshtools as mt +from pygimli.physics import ert +from pygimli.physics import traveltime as tt from pygimli.physics.petro import transFwdArchieS as ArchieTrans from pygimli.physics.petro import transFwdWyllieS as WyllieTrans -from pygimli.frameworks import PetroInversionManager, JointPetroInversionManager +from pygimli.frameworks import (PetroInversionManager, + JointPetroInversionManager) -################################################################################ +############################################################################### # We start with defining two helper functions. def createSynthModel(): @@ -89,12 +92,13 @@ def showModel(ax, model, mesh, petro=1, cMin=None, cMax=None, label=None, linewidth=0.3, color="0.2") return ax -################################################################################ + +############################################################################### # Create synthetic model # ...................... mMesh, pMesh, saturation = createSynthModel() -################################################################################ +############################################################################### # Create Petrophysical models ertTrans = ArchieTrans(rFluid=20, phi=0.3) res = ertTrans(saturation) @@ -104,7 +108,7 @@ def showModel(ax, model, mesh, petro=1, cMin=None, cMax=None, label=None, sensors = mMesh.positions()[mMesh.findNodesIdxByMarker(-99)] -################################################################################ +############################################################################### # Forward simulation # .................. # To create synthetic data sets, we assume 16 equally-spaced sensors on the @@ -112,27 +116,27 @@ def showModel(ax, model, mesh, petro=1, cMin=None, cMax=None, label=None, # complete dipole-dipole array. For the ultrasonic tomography we simulate the # travel time for every possible sensor pair. pg.info("Simulate ERT") -ERT = pg.physics.ert.ERTManager(verbose=False, sr=False) -ertScheme = pg.physics.ert.createERTData(sensors, schemeName='dd', closed=1) -ertData = ERT.simulate(mMesh, scheme=ertScheme, res=res, noiseLevel=0.01) +ertScheme = ert.createERTData(sensors, schemeName='dd', closed=1) +ertData = ert.simulate(mMesh, scheme=ertScheme, res=res, noiseLevel=0.01) pg.info("Simulate Traveltime") -TT = pg.physics.traveltime.TravelTimeManager(verbose=False) -ttScheme = pg.physics.traveltime.createRAData(sensors) -ttData = TT.simulate(mMesh, scheme=ttScheme, vel=vel, +ttScheme = tt.createRAData(sensors) +ttData = tt.simulate(mMesh, scheme=ttScheme, vel=vel, noiseLevel=0.01, noiseAbs=4e-6) -################################################################################ +############################################################################### # Conventional inversion pg.info("ERT Inversion") +ERT = ert.ERTManager(verbose=False, sr=False) resInv = ERT.invert(ertData, mesh=pMesh, zWeight=1, lam=20, verbose=False) ERT.inv.echoStatus() pg.info("Traveltime Inversion") +TT = tt.TravelTimeManager(verbose=False) velInv = TT.invert(ttData, mesh=pMesh, lam=100, useGradient=0, zWeight=1.0) TT.inv.echoStatus() -################################################################################ +############################################################################### # Petrophysical inversion (individually) pg.info("ERT Petrogeophysical Inversion") ERTPetro = PetroInversionManager(petro=ertTrans, mgr=ERT) @@ -145,7 +149,7 @@ def showModel(ax, model, mesh, petro=1, cMin=None, cMax=None, label=None, satTT = TTPetro.invert(ttData, mesh=pMesh, limits=[0., 1.], lam=5) TTPetro.inv.echoStatus() -################################################################################ +############################################################################### # Petrophysical joint inversion pg.info("Petrophysical Joint-Inversion TT-ERT") JointPetro = JointPetroInversionManager(petros=[ertTrans, ttTrans], @@ -154,7 +158,7 @@ def showModel(ax, model, mesh, petro=1, cMin=None, cMax=None, label=None, limits=[0., 1.], lam=5, verbose=False) JointPetro.inv.echoStatus() -################################################################################ +############################################################################### # Show results ERT.showData(ertData) TT.showData(ttData) diff --git a/pygimli/frameworks/methodManager.py b/pygimli/frameworks/methodManager.py index 576a9d43b..55240f230 100644 --- a/pygimli/frameworks/methodManager.py +++ b/pygimli/frameworks/methodManager.py @@ -65,7 +65,6 @@ def fit(funct, data, err=None, **kwargs): if isinstance(err, (float, int)): err = np.full(len(data), err) - mgr = ParameterInversionManager(funct, **kwargs) model = mgr.invert(data, err, **kwargs) return model, mgr.fw.response @@ -463,7 +462,6 @@ def showData(self, data=None, ax=None, **kwargs): if ax is None: fig, ax = pg.plt.subplots() - if data is None: data = self.data @@ -529,7 +527,7 @@ def showResultAndFit(self, **kwargs): """ saveFig = kwargs.pop('saveFig', None) - + if 'axs' in kwargs: axs = kwargs.pop('axs') if len(axs) != 3: @@ -538,14 +536,14 @@ def showResultAndFit(self, **kwargs): ax1 = axs[1] ax2 = axs[2] else: - fig = pg.plt.figure(figsize=kwargs.pop('figsize', (11,6))) + fig = pg.plt.figure(figsize=kwargs.pop('figsize', (11, 6))) ax = fig.add_subplot(1, 2, 1) ax1 = fig.add_subplot(2, 2, 2) ax2 = fig.add_subplot(2, 2, 4) - fig = ax.figure - + fig = ax.figure + if 'title' in kwargs: fig.suptitle(kwargs['title']) @@ -734,9 +732,17 @@ def invert(self, data=None, mesh=None, startModel=None, zWeight : float [None] Set zWeight or use defaults from regionManager. - correlationLengths - + correlationLengths : [float, float, float] + Correlation lengths for geostatistical regularization + + dip : float + rotation axis between first and last dimension (x and z) + + strike : float + rotation axis between first and second dimension (x and y) + limits: [float, float] + lower and upper value bounds for logarithmic transformation Al other are forwarded to Inversion.run From 9f30aea79d440bf0b6a2fec8de02201599b40416 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20G=C3=BCnther?= Date: Fri, 12 Aug 2022 23:38:22 +0200 Subject: [PATCH 6/8] add extractUpperSurface2dMesh --- pygimli/meshtools/__init__.py | 5 ++-- pygimli/meshtools/mesh.py | 45 +++++++++++++++++++++++++++++++++++ 2 files changed, 48 insertions(+), 2 deletions(-) diff --git a/pygimli/meshtools/__init__.py b/pygimli/meshtools/__init__.py index d07beb1e2..05d3b8ee9 100644 --- a/pygimli/meshtools/__init__.py +++ b/pygimli/meshtools/__init__.py @@ -23,12 +23,13 @@ tapeMeasureToCoordinates) from .mesh import (convert, convertMeshioMesh, convertHDF5Mesh, createMesh, - createParaMesh, createParaMesh2DGrid, createMeshFromHull, exportFenicsHDF5Mesh, exportHDF5Mesh, + createParaMesh, createParaMesh2DGrid, createMeshFromHull, + exportFenicsHDF5Mesh, exportHDF5Mesh, exportSTL, extrudeMesh, merge2Meshes, mergeMeshes, readFenicsHDF5Mesh, readGmsh, readHDF5Mesh, readHydrus2dMesh, readHydrus3dMesh, readSTL, readTetgen, readTriangle, readMeshIO, refineHex2Tet, refineQuad2Tri, - toSubsurface, fromSubsurface) + toSubsurface, fromSubsurface, extractUpperSurface2dMesh) from .polytools import createParaDomain2D # keep for backward compatibility from .polytools import (createCircle, createCube, createCylinder, createFacet, diff --git a/pygimli/meshtools/mesh.py b/pygimli/meshtools/mesh.py index 0596d2371..64a6e8be1 100644 --- a/pygimli/meshtools/mesh.py +++ b/pygimli/meshtools/mesh.py @@ -2037,5 +2037,50 @@ def createParaMesh2DGrid(sensors, paraDX=1, paraDZ=1, paraDepth=0, nLayers=11, return mesh +def extractUpperSurface2dMesh(mesh, zCut=None): + """Extract 2d mesh from the upper surface of a 3D mesh. + + Useful for showing a quick 2D plot of a 3D parameter distribution + All cell-based parameters are copied to the new mesh + + Parameters + ---------- + mesh : pg.Mesh + input mesh (3D) + zCut : float + z value to distinguish between top and bottom + + Returns + ------- + mesh2d : pg.Mesh + output 2D mesh consisting of triangles or quadrangles + + Examples + -------- + >>> import pygimli as pg + >>> from pygimli.meshtools import extractUpperSurface2dMesh + >>> mesh3d = pg.createGrid(5, 4, 3) + >>> mesh3d["val"] = pg.utils.grange(0, mesh3d.cellCount(), 1) + >>> mesh2d = extractUpperSurface2dMesh(mesh3d) + >>> ax, _ = pg.show(mesh2d, "val") + """ + if mesh.boundaryCount() == 0: + mesh.createNeighborInfos() + + zCut = zCut or pg.mean(pg.z(mesh)) + bind = [b.id() for b in mesh.boundaries() if b.outside() and + b.center().z() > zCut and b.shape().norm().z() != 0] + bMesh = mesh.createSubMesh(mesh.boundaries(bind)) + cind = np.array([mesh.boundary(i).leftCell().id() for i in bind]) + mesh2d = pg.Mesh(2) + [mesh2d.createNode(n.pos()) for n in bMesh.nodes()] + for b in bMesh.boundaries(): + mesh2d.createCell([n.id() for n in b.nodes()]) + for k in mesh.dataKeys(): + mesh2d[k] = mesh[k][cind] + + return mesh2d + + if __name__ == "__main__": pass From b3eed19e84d1ca989d279f1881041d3e57b47373 Mon Sep 17 00:00:00 2001 From: carsten-forty2 Date: Wed, 17 Aug 2022 16:28:35 +0200 Subject: [PATCH 7/8] fix JointPetro example --- pygimli/frameworks/methodManager.py | 54 ++++++++++++++++------------- pygimli/frameworks/modelling.py | 29 +++++++++------- pygimli/testing/test_show.py | 7 +++- 3 files changed, 52 insertions(+), 38 deletions(-) diff --git a/pygimli/frameworks/methodManager.py b/pygimli/frameworks/methodManager.py index 55240f230..e811e2bc3 100644 --- a/pygimli/frameworks/methodManager.py +++ b/pygimli/frameworks/methodManager.py @@ -91,13 +91,13 @@ class MethodManager(object): Attributes ---------- - verbose : bool + verbose: bool Give verbose output. - debug : bool + debug: bool Give debug output. - fop : :py:mod:`pygimli.frameworks.Modelling` + fop: :py:mod:`pygimli.frameworks.Modelling` Forward Operator instance .. knows the physics. fop is initialized by :py:mod:`pygimli.manager.MethodManager.initForwardOperator` @@ -105,7 +105,7 @@ class MethodManager(object): :py:mod:`pygimli.manager.MethodManager.createForwardOperator` method in any derived classes. - inv : :py:mod:`pygimli.frameworks.Inversion`. + inv: :py:mod:`pygimli.frameworks.Inversion`. Inversion framework instance .. knows the reconstruction approach. The attribute inv is initialized by default but can be changed overwriting @@ -273,18 +273,18 @@ def estimateError(self, data, errLevel=0.01, absError=None): Parameters ---------- - data : iterable + data: iterable Data values for which the errors should be estimated. - errLevel : float (0.01) + errLevel: float (0.01) Error level in percent/100 (i.e., 3% = 0.03). - absError : float (None) + absError: float (None) Absolute error in the unit of the data. Returns ------- - err : array + err: array Returning array of size len(data) """ if absError is not None: @@ -388,10 +388,10 @@ def invert(self, data=None, err=None, **kwargs): Parameters ---------- - dataVals : iterable + dataVals: iterable Data values to be inverted. - errVals : iterable | float + errVals: iterable | float Error value for the given data. If errVals is float we assume this means to be a global relative error and force self.estimateError to be called. @@ -421,10 +421,10 @@ def showModel(self, model, ax=None, **kwargs): Parameters ---------- - ax : mpl axes + ax: mpl axes Axes object to draw into. Create a new if its not given. - model : iterable + model: iterable Model data to be draw. Returns @@ -448,10 +448,10 @@ def showData(self, data=None, ax=None, **kwargs): Parameters ---------- - ax : mpl axes + ax: mpl axes Axes object to draw into. Create a new if its not given. - data : iterable | pg.DataContainer + data: iterable | pg.DataContainer Data values to be draw. Returns @@ -476,10 +476,10 @@ def showResult(self, model=None, ax=None, **kwargs): Parameters ---------- - ax : mpl axes + ax: mpl axes Axes object to draw into. Create a new if its not given. - model : iterable [None] + model: iterable [None] Model values to be draw. Default is self.model from the last run Returns @@ -678,7 +678,7 @@ def __init__(self, **kwargs): Attribute --------- - mesh: pg.Mesh + mesh: :gimliapi:`GIMLI::Mesh` Copy of the main mesh to be distributed to inversion and the fop. You can overwrite it with invert(mesh=mesh). """ @@ -702,12 +702,13 @@ def createMesh(self, data=None, **kwargs): def setMesh(self, mesh, **kwargs): """Set a mesh and distribute it to the forward operator""" + # keep a copy of the original mesh self.mesh = mesh self.applyMesh(mesh, **kwargs) def applyMesh(self, mesh, ignoreRegionManager=False, **kwargs): """ """ - if ignoreRegionManager: + if ignoreRegionManager is True: mesh = self.fop.createRefinedFwdMesh(mesh, **kwargs) self.fop.setMesh(mesh, ignoreRegionManager=ignoreRegionManager) @@ -718,12 +719,11 @@ def invert(self, data=None, mesh=None, startModel=None, Parameters ---------- - data : pg.DataContainer + data: pg.DataContainer - mesh : pg.Mesh [None] - - - startModel : float | iterable [None] + mesh: :gimliapi:`GIMLI::Mesh` [None] + + startModel: float | iterable [None] If set to None fop.createDefaultStartModel(dataValues) is called. @@ -882,9 +882,13 @@ def __init__(self, petro, mgr=None, **kwargs): super().__init__(fop=petrofop, **kwargs) -# Really necessary? Should a combination of petro and joint do the same class JointPetroInversionManager(MeshMethodManager): - """Joint inversion targeting at the same parameter through petrophysics.""" + """Joint inversion targeting at the same parameter through petrophysics. + + This is just syntactic sugar for the combination of + :py:mod:`pygimli.frameworks.PetroModelling` and + :py:mod:`pygimli.frameworks.JointModelling`. + """ def __init__(self, petros, mgrs): """Initialize with lists of managers and transformations""" diff --git a/pygimli/frameworks/modelling.py b/pygimli/frameworks/modelling.py index 0e479769c..afce598a3 100644 --- a/pygimli/frameworks/modelling.py +++ b/pygimli/frameworks/modelling.py @@ -609,12 +609,18 @@ def paraModel(self, model): return mod def ensureContent(self): - """""" - # Be sure the mesh is initialized when needed + """Internal function to ensure there is a valid initialized mesh. + + Initialization means the cell marker are recounted and/or there was a mesh refinement or boundary enlargement, all to fit the needs for the method depending forward problem. + """ + ## We need to call this once to be sure the mesh is initialized when needed self.mesh() def setMeshPost(self, data): - """Called when the mesh has been set successfully.""" + """Interface to be called when the mesh has been set successfully. + + Might be overwritten by child classes. + """ pass def createRefinedFwdMesh(self, mesh): @@ -654,10 +660,8 @@ def createFwdMesh_(self): self._regionChanged = False super().setMesh(m, ignoreRegionManager=True) - # This might be confusing for users: self.mesh vs. self.mesh() - # pyflakes complains about redefinition def mesh(self): - """Return mesh.""" + """Returns the current used mesh.""" self._applyRegionProperties() if self._regionManagerInUse and self._regionChanged is True: @@ -667,11 +671,14 @@ def mesh(self): def setMesh(self, mesh, ignoreRegionManager=False): """Set mesh and specify whether region manager can be ignored.""" + ### keep a copy, just in case self._baseMesh = mesh + if ignoreRegionManager is False: self._regionManagerInUse = True - if ignoreRegionManager or not self._regionManagerInUse: + ### Modelling without region manager + if ignoreRegionManager is True or not self._regionManagerInUse: self._regionManagerInUse = False if self.fop is not None: pg.critical('in use?') @@ -682,9 +689,9 @@ def setMesh(self, mesh, ignoreRegionManager=False): self.setMeshPost(mesh) return - self.clearRegionProperties() # copy the mesh to the region manager who renumber cell markers + self.clearRegionProperties() self.regionManager().setMesh(mesh) self.setDefaultBackground() @@ -804,7 +811,7 @@ def createJacobian(self, model): # pg._r("create Jacobian", self, self._jac) self.setJacobian(self._jac) # to be sure .. test if necessary - +#220817 to be changed !! # class JointModelling(Modelling): class JointModelling(MeshModelling): """Cumulative (joint) forward operator.""" @@ -857,9 +864,7 @@ def setMesh(self, mesh, **kwargs): # to be removed from here for fi in self.fops: fi.setMesh(mesh) - self.setRegionManager(self.fops[0].regionManagerRef()) - - +#220817 to be implemented!! # class JointMeshModelling(JointModelling): # def __init__(self, fopList): # super().__init__(self, fopList) diff --git a/pygimli/testing/test_show.py b/pygimli/testing/test_show.py index c7f9b9b4d..e04c66255 100644 --- a/pygimli/testing/test_show.py +++ b/pygimli/testing/test_show.py @@ -208,9 +208,14 @@ def testCBarLevels(): def testShowPV(): """ - # Run from notebook import pygimli as pg from pygimli.testing.test_show import testShowPV + + import pyvista as pv + print(pv.__version__) + import panel as pnl + print(pnl.__version__) + testShowPV() """ m1 = mt.createCube() From 7c47375aab71e0fc2d03f55d5cef062bc1c6c412 Mon Sep 17 00:00:00 2001 From: carsten-forty2 Date: Thu, 18 Aug 2022 18:58:22 +0200 Subject: [PATCH 8/8] last fixes before 1.3 --- core/python/generate_pygimli_code.py | 9 +- core/src/gimli.h | 75 +++++++++----- core/src/mesh.cpp | 141 +++++++++++++++++++++------ core/src/meshentities.cpp | 11 ++- core/src/meshentities.h | 4 + core/src/node.cpp | 9 ++ core/src/node.h | 4 +- core/src/plane.cpp | 20 +++- core/src/plane.h | 14 ++- core/src/vector.cpp | 12 +++ core/src/vector.h | 4 + pygimli/meshtools/mesh.py | 13 ++- pygimli/meshtools/polytools.py | 24 +++-- pygimli/testing/test_polytools.py | 78 ++++++++++++++- pygimli/viewer/pv/utils.py | 10 ++ 15 files changed, 343 insertions(+), 85 deletions(-) diff --git a/core/python/generate_pygimli_code.py b/core/python/generate_pygimli_code.py index f63b93835..c13d61f98 100644 --- a/core/python/generate_pygimli_code.py +++ b/core/python/generate_pygimli_code.py @@ -373,7 +373,11 @@ def generate(defined_symbols, extraIncludes): 'getNonEmptyRow', 'getSubstrings', 'abs', - 'type'] + 'type', + '::GIMLI::print', + 'print', + 'range', + ] ) exclude(main_ns.free_operators, @@ -455,6 +459,9 @@ def generate(defined_symbols, extraIncludes): 'IVectorIter', 'RVectorIter', 'FunctionDD', + '::GIMLI::print', + 'print', + 'range', ] for c in main_ns.free_functions(): diff --git a/core/src/gimli.h b/core/src/gimli.h index c4022a444..82119c76f 100644 --- a/core/src/gimli.h +++ b/core/src/gimli.h @@ -126,12 +126,62 @@ typedef int64_t int64; #define __FILENAME__ __FILE__ #endif + +#ifndef PYGIMLI_CAST // castxml complains on older gcc/clang +inline std::string str(){ return "";} +#endif +//! General template for conversion to string, should supersede all sprintf etc. +template< typename T > inline std::string str(const T & v){ + std::ostringstream os; + os << v; + return os.str(); +} +enum LogType {Verbose, Info, Warning, Error, Debug, Critical}; +DLLEXPORT void log(LogType type, const std::string & msg); + + +template +std::string str(Value v, Values... vs){ + std::ostringstream os; + using expander = int[]; + os << v; // first + (void) expander{ 0, (os << " " << vs, void(), 0)... }; + return os.str(); +} +template +void log(LogType type, Values... vs){ + return log(type, str(vs...)); +} + +// simple python like std out +inline void print() { + std::cout << std::endl; +} + +template +void print(std::ostream& s, Head&& head) { + s << std::forward(head) << std::endl; +} + +template +void print(std::ostream& s, Head&& head, Tail&&... tail) { + s << std::forward(head) << " "; + print(s, std::forward(tail)...); +} + +template +void print(Args&&... args) { + print(std::cout, std::forward(args)...); +} + #define WHERE GIMLI::str(__FILENAME__) + ":" + GIMLI::str(__LINE__) + "\t" #define WHERE_AM_I WHERE + "\t" + GIMLI::str(__ASSERT_FUNCTION) + " " #define TO_IMPL WHERE_AM_I + " not yet implemented\n " + GIMLI::versionStr() + "\nPlease send the messages above, the commandline and all necessary data to the author." #define __M std::cout << "*** " << WHERE << std::endl; #define __MS(str) std::cout << "*** " < #endif -inline std::string str(){ return "";} -//! General template for conversion to string, should supersede all sprintf etc. -template< typename T > inline std::string str(const T & v){ - std::ostringstream os; - os << v; - return os.str(); -} -enum LogType {Verbose, Info, Warning, Error, Debug, Critical}; -DLLEXPORT void log(LogType type, const std::string & msg); - -template -std::string str(Value v, Values... vs){ - std::ostringstream os; - using expander = int[]; - os << v; // first - (void) expander{ 0, (os << " " << vs, void(), 0)... }; - return os.str(); -} -#ifndef PYGIMLI_CAST // castxml complains on older gcc/clang -#endif -template -void log(LogType type, Values... vs){ - return log(type, str(vs...)); -} - DLLEXPORT std::string versionStr(); DLLEXPORT std::string authors(); diff --git a/core/src/mesh.cpp b/core/src/mesh.cpp index 22cbd6b83..809be2284 100644 --- a/core/src/mesh.cpp +++ b/core/src/mesh.cpp @@ -19,12 +19,12 @@ #include "mesh.h" #include "kdtreeWrapper.h" +#include "line.h" #include "memwatch.h" #include "meshentities.h" #include "node.h" -#include "line.h" +#include "plane.h" #include "shape.h" - #include "sparsematrix.h" #include "stopwatch.h" @@ -178,6 +178,10 @@ Node * Mesh::createNodeGC_(const RVector3 & pos, int marker){ Index oldCount = this->nodeCount(); Node *n = this->createNodeWithCheck(pos); n->setMarker(marker); + + if ((this->nodeCount() == oldCount)){ + if (n->state() == No) n->setState(Original); + } if ((this->dim() == 3) and (this->nodeCount() > oldCount)){ @@ -501,15 +505,18 @@ Boundary * findSecParent(const std::vector < Node * > & v){ } return 0; } + Boundary * Mesh::copyBoundary(const Boundary & bound, double tol, bool check){ + bool debug = false; + #define _D(...) if (debug) __MSP(__VA_ARGS__) + _D("copyBoundary", bound.ids()) std::vector < Node * > nodes(bound.nodeCount()); - bool isFreeFace = false; - // bool isSubFace = false; - + bool isFreeFace = false; //** the new face is no subface + + std::vector < Node * > oldNodes; std::vector < Node * > conNodes; std::vector < Node * > secNodes; - std::vector < Node * > subNodes; // __M for (Index i = 0; i < nodes.size(); i ++) { @@ -519,87 +526,159 @@ Boundary * Mesh::copyBoundary(const Boundary & bound, double tol, bool check){ nodes[i] = createNode(bound.node(i).pos(), tol); } else { // 3D poly tests fail!! .. need to be checked and fixed TODO - nodes[i] = createNodeWithCheck(bound.node(i).pos(), tol); + if (check== true){ + nodes[i] = createNodeWithCheck(bound.node(i).pos(), tol); + } else { + nodes[i] = createNode(bound.node(i).pos()); + } } nodes[i]->setMarker(bound.node(i).marker()); // __MS(nodes[i]->state()) switch (nodes[i]->state()){ case NodeState::No: - // at least one node is not in boundary + // at least one node is not inside boundary isFreeFace = true; break; case NodeState::Secondary: secNodes.push_back(nodes[i]); break; // __MS(*nodes[i]) case NodeState::Connected: conNodes.push_back(nodes[i]); break; + case NodeState::Original: + oldNodes.push_back(nodes[i]); + break; + } } - Boundary * b = 0; Boundary * parent = 0; + Boundary * ret = 0; - if (bound.rtti() == MESH_POLYGON_FACE_RTTI){ + std::vector < Node * > subNodes; // nodes for the new subface + if (bound.rtti() == MESH_POLYGON_FACE_RTTI && check == true){ + + _D("connectedNodes:", conNodes, + "secondaryNodes:", secNodes, + "origNodes:", oldNodes) + _D("new face nodes:", nodes, "is subface", not isFreeFace) - Boundary * conParent = findBoundary(conNodes); Boundary * secParent = findSecParent(secNodes); + //** identify if face is freeface if (!isFreeFace){ - conParent = findBoundary(conNodes); - + + std::set < Boundary * > conParentCand = findBoundaries(conNodes); + Boundary * conParent = 0; + if (conParentCand.size() > 0 ){ + for (auto *b: conParentCand){ + if (b->shape().plane().compare(bound.shape().plane(), + TOLERANCE, true)){ + conParent = b; + break; + } + } + } if (conNodes.size() && secNodes.size()){ - if (conParent != secParent){ + + if (!conParent || conParent != secParent){ isFreeFace = true; } - subNodes = conNodes; - subNodes.insert(subNodes.end(), secNodes.begin(), secNodes.end()); + if (conParent){ + _D("conParent", conParent->ids()) + } else { + _D("no parent for connected nodes") + } + if (secParent){ + _D("secParent", secParent->ids()) + } else { + _D("no parent for secondary nodes") + } + + // subNodes = conNodes; + // subNodes.insert(subNodes.end(), secNodes.begin(), secNodes.end()); + subNodes = nodes; + _D("subNodes", subNodes) parent = secParent; - } - if (conNodes.size()){ + } else if (conNodes.size()){ if (!conParent){ isFreeFace = true; + _D("no conParent") + } else { + _D("conParent", conParent->ids()) + + if (oldNodes.size()){ + //original nodes may not be part of connected Parent + for (auto *n: oldNodes){ + if (!conParent->shape().touch(n->pos())) { + _D("node not on connected parent", n->id()) + isFreeFace = true; + } + } + } } - subNodes = conNodes; + subNodes = nodes; parent = conParent; - } - if (secNodes.size()){ + } else if (secNodes.size()){ if(!secParent){ isFreeFace = true; + _D("no secondary parent") + } else if (!secParent->shape().plane().compare( + bound.shape().plane(), TOLERANCE, true)){ + + _D("secParent", secParent->ids()) + _D("secParent.Plane:", secParent->shape().plane()); + _D("bound.Plane:", bound.shape().plane()); + isFreeFace = true; } subNodes = secNodes; parent = secParent; + } else if (oldNodes.size()){ + isFreeFace = true; } } - + + //** create new face + _D("is subface", not isFreeFace, subNodes.size(), nodes.size()) + if (isFreeFace){ - b = createBoundaryChecked_< PolygonFace >(nodes, - bound.marker(), check); - parent = b; + ret = createBoundaryChecked_< PolygonFace >(nodes, + bound.marker(), check); + _D("added freeFace: ", ret->ids()) + parent = ret; } else { if (subNodes.size() > 2){ if (parent){ for (auto *n: secNodes){ parent->delSecondaryNode(n); + n->setState(No); + } + auto *p = dynamic_cast< PolygonFace * >(parent); + _D("addSubface: ", p->subfaceCount()) + p->addSubface(subNodes); + _D("subfacecount: ", p->subfaceCount()) + + if (debug) for (auto i: range(p->subfaceCount())){ + _D("added subface:", p->subface(i)) } - dynamic_cast< PolygonFace * >(parent)->addSubface(subNodes); } else { log(Error, "no parent boundary"); } - b = parent; + ret = parent; + } else { + log(Error, "new subface but only two nodes"); } } const PolygonFace & f = dynamic_cast< const PolygonFace & >(bound); if (f.subfaceCount() > 0){ log(Error, "Can't yet copy a boundary with subfaces"); } - for (Index i = 0; i < f.holeMarkers().size(); i ++ ){ dynamic_cast< PolygonFace* >(parent)->addHoleMarker( f.holeMarkers()[i]); } } else { // if no Polygonface - b = createBoundary(nodes, bound.marker(), check); + ret = createBoundary(nodes, bound.marker(), check); } - return b; + return ret; } Node & Mesh::node(Index i) { @@ -1998,7 +2077,6 @@ Mesh Mesh::createSubMesh(const std::vector< Node * > & nodes) const { return mesh; } - void Mesh::addData(const std::string & name, const RVector & data) { // std::cout << "add export Data: " << name << " " << min(data) << " " << max(data) << std::endl; if (dataMap_.count(name)){ @@ -2093,6 +2171,7 @@ void Mesh::mapBoundaryMarker(const std::map < int, int > & aMap){ } } } + void Mesh::prolongateEmptyCellsValues(RVector & vals, double background) const { IndexArray emptyList(find(abs(vals) < TOLERANCE)); if (emptyList.size() == 0) return; diff --git a/core/src/meshentities.cpp b/core/src/meshentities.cpp index eb46f2c27..212dbb373 100644 --- a/core/src/meshentities.cpp +++ b/core/src/meshentities.cpp @@ -18,10 +18,10 @@ #include "meshentities.h" -#include "node.h" -#include "shape.h" #include "line.h" #include "mesh.h" +#include "node.h" +#include "shape.h" #include #include @@ -73,14 +73,17 @@ Boundary * findBoundary(const std::vector < Node * > & n) { else if (n.size() == 3) return findBoundary(*n[0], *n[1], *n[2]); else if (n.size() == 4) return findBoundary(*n[0], *n[1], *n[2], *n[3]); + return findBoundary_(findBoundaries(n)); +} + +std::set < Boundary * > findBoundaries(const std::vector < Node * > & n){ std::vector < std::set< Boundary * > > bs(n.size()); for (uint i = 0; i < n.size(); i ++) bs[i] = n[i]->boundSet(); std::set < Boundary * > common; intersectionSet(common, bs); - - return findBoundary_(common); + return common; } Boundary * findCommonBoundary(const Cell & c1, const Cell & c2){ diff --git a/core/src/meshentities.h b/core/src/meshentities.h index 076e55e6c..1aaefe37c 100644 --- a/core/src/meshentities.h +++ b/core/src/meshentities.h @@ -43,6 +43,10 @@ DLLEXPORT Boundary * findBoundary(const Node & n1, const Node & n2, const Node & DLLEXPORT Boundary * findBoundary(const std::vector < Node * > & n); +/*! Find all boundaries that have these nodes in common. */ +DLLEXPORT std::set < Boundary * > +findBoundaries(const std::vector < Node * > & n); + /*! */ DLLEXPORT Boundary * findCommonBoundary(const Cell & c1, const Cell & c2); diff --git a/core/src/node.cpp b/core/src/node.cpp index 02eff45b3..817a06f9e 100644 --- a/core/src/node.cpp +++ b/core/src/node.cpp @@ -30,6 +30,15 @@ std::ostream & operator << (std::ostream & str, const GIMLI::Node & n){ return str; } +std::ostream & operator << (std::ostream & str, + const std::vector < GIMLI::Node * > & nodes){ + for (auto *n: nodes) { + str << n->id() << " "; + } + str << std::endl; + return str; +} + Node::Node(){ init_(); marker_ = 0; diff --git a/core/src/node.h b/core/src/node.h index ee52bf302..aa983b8d1 100644 --- a/core/src/node.h +++ b/core/src/node.h @@ -28,7 +28,7 @@ namespace GIMLI{ -enum NodeState{No, Secondary, Connected}; +enum NodeState{No, Original, Secondary, Connected}; //! 3D Node /*! @@ -148,6 +148,8 @@ class DLLEXPORT Node : public BaseEntity { DLLEXPORT std::ostream & operator << (std::ostream & str, const GIMLI::Node & node); +DLLEXPORT std::ostream & operator << (std::ostream & str, const std::vector < GIMLI::Node * > & nodes); + inline bool operator == (const Node & n1, const Node & n2) { return n1.pos() == n2.pos(); } } // namespace GIMLI diff --git a/core/src/plane.cpp b/core/src/plane.cpp index 46e158b7f..c44d95dd2 100644 --- a/core/src/plane.cpp +++ b/core/src/plane.cpp @@ -187,11 +187,23 @@ Line Plane::intersect(const Plane & plane, double tol){ return Line(lineP0, lineP1); } -bool Plane::compare(const Plane & plane, double tol){ - if (norm_.distance(plane.norm()) < tol && - std::fabs(d_ - plane.d()) < tol) return true; +bool Plane::compare(const Plane & plane, double tol, bool bothDirs){ + if (bothDirs == true){ + // __MS(std::fabs(std::fabs(d_) - std::fabs(plane.d()))) + // __MS(norm_.distance(plane.norm())) + // __MS(std::fabs(norm_.distance(plane.norm()) - 2.0)) + + if (std::fabs(std::fabs(d_) - std::fabs(plane.d())) > tol) return false; + + if (norm_.distance(plane.norm()) < tol || + std::fabs(norm_.distance(plane.norm()) - 2.0) < tol) return true; + + } else { + if (std::fabs(d_ - plane.d()) > tol) return false; + if (norm_.distance(plane.norm()) < tol) return true; + } + return false; } - } //namespace GIMLI diff --git a/core/src/plane.h b/core/src/plane.h index fc4489ba8..91d238c24 100644 --- a/core/src/plane.h +++ b/core/src/plane.h @@ -54,14 +54,20 @@ class DLLEXPORT Plane{ Plane & operator = (const Plane & plane); /*! Equal_to operator */ - inline bool operator == (const Plane & plane){ return this->compare(plane); } + inline bool operator == (const Plane & plane){ + return this->compare(plane); + } /*! Not_equal_to operator */ - inline bool operator != (const Plane & plane){ return !(*this== plane); } + inline bool operator != (const Plane & plane){ + return !(*this== plane); + } /*! Compare two planes with a given tolerance. Check if both norms and distances are equal. - | norm - p.norm | < tol && | d_ - p.d | < tol */ - bool compare(const Plane & p, double tol = TOLERANCE); + | norm - p.norm | < tol && | d_ - p.d | < tol. + If bothDirection is True | norm - p.norm | can be 2.0 to allow for coplanar check with different orientations.*/ + bool compare(const Plane & p, double tol=TOLERANCE, bool bothDirs=false); + /*! Returns true if the plane is valid and pos touch this plane. Touch when plane.distance(pos) < tol. */ diff --git a/core/src/vector.cpp b/core/src/vector.cpp index 73fd92ff4..10b2dd366 100644 --- a/core/src/vector.cpp +++ b/core/src/vector.cpp @@ -93,6 +93,18 @@ void Vector< RVector3 >::clean(){ } } +IndexArray range(Index start, Index stop, Index step){ + IndexArray ret(0); + for (Index i = start; i < stop; i += step){ + ret.push_back(i); + } + return ret; +} + +IndexArray range(Index stop){ + return range(0, stop, 1); +} + } // namespace GIMLI{ diff --git a/core/src/vector.h b/core/src/vector.h index a6797dd56..2158ef80a 100644 --- a/core/src/vector.h +++ b/core/src/vector.h @@ -64,6 +64,10 @@ template < class ValueType, class A > class __VectorExpr; IndexArray find(const BVector & v); +DLLEXPORT IndexArray range(Index start, Index stop, Index step=1); +DLLEXPORT IndexArray range(Index stop); + + #ifndef PYGIMLI_CAST inline void Dump(const void * mem, unsigned int n) { const char * p = reinterpret_cast< const char *>(mem); diff --git a/pygimli/meshtools/mesh.py b/pygimli/meshtools/mesh.py index 64a6e8be1..42974bef2 100644 --- a/pygimli/meshtools/mesh.py +++ b/pygimli/meshtools/mesh.py @@ -2045,14 +2045,14 @@ def extractUpperSurface2dMesh(mesh, zCut=None): Parameters ---------- - mesh : pg.Mesh - input mesh (3D) - zCut : float + mesh: :gimliapi:`GIMLI::Mesh` + Input mesh (3D) + zCut: float z value to distinguish between top and bottom Returns ------- - mesh2d : pg.Mesh + mesh2d: :gimliapi:`GIMLI::Mesh` output 2D mesh consisting of triangles or quadrangles Examples @@ -2071,11 +2071,14 @@ def extractUpperSurface2dMesh(mesh, zCut=None): bind = [b.id() for b in mesh.boundaries() if b.outside() and b.center().z() > zCut and b.shape().norm().z() != 0] bMesh = mesh.createSubMesh(mesh.boundaries(bind)) - cind = np.array([mesh.boundary(i).leftCell().id() for i in bind]) + mesh2d = pg.Mesh(2) [mesh2d.createNode(n.pos()) for n in bMesh.nodes()] for b in bMesh.boundaries(): mesh2d.createCell([n.id() for n in b.nodes()]) + + # copy data + cind = np.array([mesh.boundary(i).leftCell().id() for i in bind]) for k in mesh.dataKeys(): mesh2d[k] = mesh[k][cind] diff --git a/pygimli/meshtools/polytools.py b/pygimli/meshtools/polytools.py index 20a9b4f67..d15081a10 100644 --- a/pygimli/meshtools/polytools.py +++ b/pygimli/meshtools/polytools.py @@ -1081,7 +1081,8 @@ def createParaMeshPLC(sensors, paraDX=1, paraDepth=-1, paraBoundary=2, def createParaMeshSurface(sensors, paraBoundary=None, boundary=-1, - surfaceMeshQuality=30, addTopo=None): + surfaceMeshQuality=30, surfaceMeshArea=0, + addTopo=None): r"""Create surface mesh for an 3D inversion parameter mesh. Topographic information (non-zero z-coodinate) can be from sensors @@ -1100,9 +1101,12 @@ def createParaMeshSurface(sensors, paraBoundary=None, boundary=-1, Boundary width to be appended for domain prolongation in relative para domain size. - surfaceMeshQuality: float (30) + surfaceMeshQuality: float [30] Quality of the surface mesh. + surfaceMeshArea: float [0] + Max cell Size for parametric domain. + addTopo: [[x,y,z],] Number of additional nodes for topography. @@ -1164,7 +1168,9 @@ def createParaMeshSurface(sensors, paraBoundary=None, boundary=-1, # find parameter extent paraRect = pg.meshtools.createRectangle(pnts=sensors[:, 0:2], minBB=True, - minBBOffset=paraBoundary) + minBBOffset=paraBoundary, + area=surfaceMeshArea, + ) for i in range(4): paraRect.boundary(i).setMarker(i + 5) paraRect.node(i).setMarker((i % 4 + 5) * 10) @@ -1188,7 +1194,8 @@ def createParaMeshSurface(sensors, paraBoundary=None, boundary=-1, def createParaMeshPLC3D(sensors, paraDX=0, paraDepth=-1, paraBoundary=None, paraMaxCellSize=0.0, boundary=None, - boundaryMaxCellSize=0, surfaceMeshQuality=30, + boundaryMaxCellSize=0, + surfaceMeshQuality=30, surfaceMeshArea=0, addTopo=None, isClosed=False, **kwargs): r"""Create a geometry (PLC) for an 3D inversion parameter mesh. @@ -1224,8 +1231,11 @@ def createParaMeshPLC3D(sensors, paraDX=0, paraDepth=-1, paraBoundary=None, Boundary width to be appended for domain prolongation in relative para domain size. - surfaceMeshQuality: float (30) + surfaceMeshQuality: float [30] Quality of the surface mesh. + + surfaceMeshArea: float [0] + Max boundary size for surface area in parametric region. addTopo: [[x,y,z],] Number of additional nodes for topography. @@ -1245,7 +1255,9 @@ def createParaMeshPLC3D(sensors, paraDX=0, paraDepth=-1, paraBoundary=None, surface = pg.meshtools.createParaMeshSurface( sensors, paraBoundary=paraBoundary, boundary=boundary, - surfaceMeshQuality=surfaceMeshQuality, addTopo=addTopo) + surfaceMeshQuality=surfaceMeshQuality, + surfaceMeshArea=surfaceMeshArea, + addTopo=addTopo) # find depth and paradepth xSpan = (max(sensors[:, 0]) - min(sensors[:, 0])) diff --git a/pygimli/testing/test_polytools.py b/pygimli/testing/test_polytools.py index f99fcc51b..0f836d1d0 100644 --- a/pygimli/testing/test_polytools.py +++ b/pygimli/testing/test_polytools.py @@ -300,8 +300,10 @@ def test_smallcube_in_bigcube(self): def test_cyl_on_cyl(self): # merge only works if smaller face merged into larger face on contact plane segs = 12 - c1 = mt.createCylinder(radius=2, nSegments=segs, boundaryMarker=1) - c2 = mt.createCylinder(radius=1, nSegments=segs, boundaryMarker=2) + c1 = mt.createCylinder(radius=2, marker=1, + nSegments=segs, boundaryMarker=1) + c2 = mt.createCylinder(radius=1, marker=2, + nSegments=segs, boundaryMarker=2) c1.translate([0, 0, 0.5]) c2.translate([0, 0, -0.5]) @@ -311,9 +313,12 @@ def test_cyl_on_cyl(self): self.assertEqual(w.boundaryCount(), segs*2 + 3) # w.exportBoundaryVTU('w') - # pg.show(w) + #pg.show(w) + m = mt.createMesh(w) + pg.show(m, m.cellMarkers()) # pg.show(mt.createMesh(w)) - + + def test_face_in_face(self): """Test subface with different marker constructed with hole marker.""" w = mt.createCube(marker=1, boundaryMarker=1) @@ -348,6 +353,70 @@ def test_face_in_face(self): # mesh.exportBoundaryVTU('b.vtu') #pg.show(mesh) + + def test_cube_cube_halfside(self): + """Add half size cube on another cube""" + D=1 + H=1 + W=3 + c1 = mt.createCube([D, D, H/2], pos=[0.0, 0.0, H/4]) + c1L = pg.Mesh(c1) + c1L.addRegionMarker([0.0, 0.0, 0.0], marker=2) + # c1R = pg.Mesh(c1L) + # c1R.translate([0.0, W-D, 0.0]) + + c2 = mt.createCube([D, W-2*D, H], pos=[0, -(W)/2+D/2-0.00, 0], marker=2) + #plc = mt.merge([c2, c1L, c1R]) + plc = mt.merge([c2, c1L]) + #plc.exportPLC('cubecut') + pg.show(plc) + + self.assertEqual(plc.nodeCount(), 14) + self.assertEqual(plc.boundaryCount(), 11) + + m = mt.createMesh(plc) + #m = mt.createMesh(plc, tetgen='tetgen-1.4.3') + pg.show(m, m.cellMarkers()) + + + def test_cube_cut_cube(self): + """Test cube cut from cube on two neighbour faces.""" + c1 = mt.createCube(marker=1) + c2 = mt.createCube(marker=2) + c2.scale([0.5, 0.5, 0.5]) + c2.translate([0.25, -0.25, 0.25]) + + plc = mt.mergePLC3D([c1, c2]) + + self.assertEqual(plc.nodeCount(), 15) + self.assertEqual(plc.boundaryCount(), 9) + + c3 = mt.createCube(isHole=True, marker=3) + c3.scale([0.3, 0.3, 0.3]) + c3.translate([0.35, 0.35, 0.35]) + plc = mt.mergePLC3D([plc, c3]) + + c4 = mt.createCube(isHole=False, marker=4) + c4.scale([0.2, 0.2, 1.0]) + c4.translate([-0.4, 0.0, 0.0]) + plc = mt.mergePLC3D([plc, c4]) + + c5 = mt.createCube(isHole=False, marker=5) + c5.scale([0.8, 0.2, 0.2]) + c5.translate([0.1, 0.0, -0.4]) + plc = mt.mergePLC3D([c1, c4, c5]) + + plc.exportPLC('cubecut') + pg.show(plc) + + #m = mt.createMesh(plc, tetgen='tetgen-1.4.3') + m = mt.createMesh(plc) + + pg.show(m, m.cellMarkers()) + + m.exportVTK('cubecut') + + def test_face_inCube(self): # plc = mt.createCube() @@ -376,5 +445,6 @@ def test_appendTetrahedron(self): # pass + if __name__ == '__main__': unittest.main() diff --git a/pygimli/viewer/pv/utils.py b/pygimli/viewer/pv/utils.py index 2b546c3f9..df905371a 100644 --- a/pygimli/viewer/pv/utils.py +++ b/pygimli/viewer/pv/utils.py @@ -38,6 +38,12 @@ def pgMesh2pvMesh(mesh, data=None, label=None, boundaries=False): Parameter to distribute to cells/nodes. """ if boundaries: + mesh.createNeighbourInfos() + + if mesh.cellCount() == 0: + ### mesh is already a boundary mesh + return pgMesh2pvMesh(mesh, data, label) + b = mesh.createSubMesh(mesh.boundaries([b.id() for b in mesh.boundaries() if b.outside() or b.marker() != 0])) return pgMesh2pvMesh(b, data, label) @@ -83,6 +89,10 @@ def pgMesh2pvMesh(mesh, data=None, label=None, boundaries=False): if label is None: label = 'Cell data' grid.cell_data[label] = np.asarray(data) + elif len(data) == mesh.boundaryCount(): + if label is None: + label = 'Boundary data' + grid.cell_data[label] = np.asarray(data) elif len(data) == mesh.nodeCount(): if label is None: label = 'Node data'