"""
Part of comet/pyhed/misc
"""
# COMET - COupled Magnetic resonance Electrical resistivity Tomography
# Copyright (C) 2019 Nico Skibbe
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
import pygimli as pg
import numpy as np
def __createH2(inmesh):
""" Creates a H2 refined mesh. """
inmesh_cp = pg.Mesh(inmesh)
outmesh = inmesh_cp.createH2()
outmesh.createNeighbourInfos()
return outmesh
def __intergrationMat(mesh1, mesh2, order=1):
""" Creates an integration matrix for H2 refined meshs. """
if mesh1.cellCount() > mesh2.cellCount():
coarse = mesh2
fine = mesh1
else:
coarse = mesh1
fine = mesh2
mat = pg.core.SparseMapMatrix(r=coarse.cellCount(), c=fine.cellCount())
for cell in fine.cells():
cid = coarse.findCell(cell.center()).id()
mat.addVal(cid, cell.id(), 1./4**order)
return mat
def __intergrationMatMixed(mesh1, mesh2):
""" Creates an integration matrix for H2 refined meshs. """
if mesh1.cellCount() > mesh2.cellCount():
coarse = mesh2
fine = mesh1
else:
coarse = mesh1
fine = mesh2
mat = pg.core.SparseMapMatrix(r=coarse.cellCount(), c=fine.cellCount())
cids = {}
for cell in fine.cells():
# for each fine cell find corresponding coarse cell
coarse_id = coarse.findCell(cell.center()).id()
cidsofar = cids.get(coarse_id, [])
cidsofar.append(cell.id())
# save ids of fine cell for each coarse cell
cids[coarse_id] = cidsofar
for ccid in range(coarse.cellCount()):
# get all ids of fine cells in coarse cell
fine_ids = cids[ccid]
# how many?
factor = len(fine_ids)
for fcid in fine_ids:
# add entry in integration mat, that will sum up to 1 in the end
# 1 is original cell, 4 for H1, 16 for H2...
mat.addVal(ccid, fcid, 1./factor)
return mat
[docs]def createH2(inmesh, order=1, integration_mat=False):
""" Creates a H-2**N refined mesh with order N and Integration matrix. """
mesh = pg.Mesh(inmesh)
for i in range(order):
mesh = __createH2(mesh)
if integration_mat:
mat = __intergrationMat(inmesh, mesh, order=order)
return (mesh, mat)
else:
return mesh
[docs]def sameGeometry(mesh1, mesh2, atol=1e-8, rtol=1e-5):
# first some short cuts to safe calculation time
# same memory: reference
if mesh1 is mesh2:
return True
if mesh1.nodeCount() != mesh2.nodeCount():
return False
if mesh1.cellCount() != mesh2.cellCount():
return False
# extensive testing
if not np.allclose(np.array(mesh1.positions()),
np.array(mesh2.positions()),
atol=atol, rtol=rtol):
return False
cells1 = np.array([np.array(cell1.ids()) for cell1 in mesh1.cells()])
cells2 = np.array([np.array(cell2.ids()) for cell2 in mesh2.cells()])
if not np.array_equiv(cells1, cells2):
return False
# further assuming same nodes and cells results in same boundaries
# and neglecting markers:
return True
[docs]def createConstraintMesh(mesh):
"""
Creates a refined mesh, where every triangle is divided in three
triangles, with the midpoint of the original cells as new node.
This is only done for each boundary that has a right and a left cell
(no boundary edges) and is usefull for displaing the boundary constraints
of the original mesh.
Example:
--------
>>> import pygimli as pg
>>> import numpy as np
>>> mesh = pg.load('invmesh.bms')
>>> cw = np.load('constraints.npy')
>>> cmesh = createConstraintMesh(mesh)
>>> pg.show(cmesh, data=cw[np.array(cmesh.cellMarkers(), dtype=int)])
"""
# need left, right cell info
mesh.createNeighbourInfos()
# output mesh
cmesh = pg.Mesh(mesh.dim())
# do for all bounds: build left cell to midpoint, same for right if any
i = 0
for bound in mesh.boundaries():
# right, if not mesh boundary
right_cell = bound.rightCell()
if right_cell is not None:
n0 = cmesh.createNodeWithCheck(bound.node(0).pos())
n1 = cmesh.createNodeWithCheck(bound.node(1).pos())
# only then there is a constraint defined
# left
n2 = cmesh.createNodeWithCheck(bound.leftCell().center())
# right
n3 = cmesh.createNodeWithCheck(right_cell.center())
cmesh.createTriangle(n0, n1, n2, i)
cmesh.createTriangle(n0, n1, n3, i)
i += 1
# finish mesh
cmesh.createNeighbourInfos()
return cmesh
# The End