simpeg/discretize

Discretization tools for finite volume and inverse problems.

SimPEG

Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications.

The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features:

• modular with respect to the spacial discretization, optimization routine, and geophysical problem
• built with the inverse problem in mind
• provides a framework for geophysical and hydrogeologic problems
• supports 1D, 2D and 3D problems
• designed for large-scale inversions

You are welcome to join our forum and engage with people who use and develop SimPEG at: https://simpeg.discourse.group/.

Weekly meetings are open to all. They are generally held on Wednesdays at 10:30am PDT. Please see the calendar (GCAL, ICAL) for information on the next meeting.

Overview Video

Working towards all the Geophysics, but Backwards - SciPy 2016

Citing SimPEG

There is a paper about SimPEG!

Cockett, R., Kang, S., Heagy, L. J., Pidlisecky, A., & Oldenburg, D. W. (2015). SimPEG: An open source framework for simulation and gradient based parameter estimation in geophysical applications. Computers & Geosciences.

BibTex:

@article{cockett2015simpeg,
  title={SimPEG: An open source framework for simulation and gradient based parameter estimation in geophysical applications},
  author={Cockett, Rowan and Kang, Seogi and Heagy, Lindsey J and Pidlisecky, Adam and Oldenburg, Douglas W},
  journal={Computers \& Geosciences},
  year={2015},
  publisher={Elsevier}
}

Electromagnetics

If you are using the electromagnetics module of SimPEG, please cite:

Lindsey J. Heagy, Rowan Cockett, Seogi Kang, Gudni K. Rosenkjaer, Douglas W. Oldenburg, A framework for simulation and inversion in electromagnetics, Computers & Geosciences, Volume 107, 2017, Pages 1-19, ISSN 0098-3004, http://dx.doi.org/10.1016/j.cageo.2017.06.018.

BibTex:

@article{heagy2017,
    title= "A framework for simulation and inversion in electromagnetics",
    author= "Lindsey J. Heagy and Rowan Cockett and Seogi Kang and Gudni K. Rosenkjaer and Douglas W. Oldenburg",
    journal= "Computers & Geosciences",
    volume = "107",
    pages = "1 - 19",
    year = "2017",
    note = "",
    issn = "0098-3004",
    doi = "http://dx.doi.org/10.1016/j.cageo.2017.06.018"
}

Questions

If you have a question regarding a specific use of SimPEG, the fastest way to get a response is by posting on our Discourse discussion forum: https://simpeg.discourse.group/. Alternatively, if you prefer real-time chat, you can join our Mattermost Team at https://mattermost.softwareunderground.org/simpeg. Please do not create an issue to ask a question.

Meetings

SimPEG hosts weekly meetings for users to interact with each other, for developers to discuss upcoming changes to the code base, and for discussing topics related to geophysics in general. Currently our meetings are held every Wednesday, alternating between a mornings (10:30 am pacific time) and afternoons (3:00 pm pacific time) on even numbered Wednesdays. Find more info on our Mattermost.

Links

Website: https://simpeg.xyz

Forums: https://simpeg.discourse.group/

Mattermost (real time chat): https://mattermost.softwareunderground.org/simpeg

Documentation: https://docs.simpeg.xyz

Code: https://github.com/simpeg/simpeg

Tests: https://dev.azure.com/simpeg/simpeg/_build

Bugs & Issues: https://github.com/simpeg/simpeg/issues

Contributing

We always welcome contributions towards SimPEG whether they are adding new code, suggesting improvements to existing codes, identifying bugs, providing examples, or anything that will improve SimPEG. Please checkout the contributing guide for more information on how to contribute.

discretize's Issues

Null Space of curl curl operation

From @rowanc1 on November 4, 2015 18:49

From @lheagy on June 29, 2015 21:24

When frequency is small, there is a null space that can cause issues. It can be corrected by enforcing conditions on the divergence of the field and or flux at play (not that if solving the 2D cyl mesh problem for h or e, this isn't an issue)

Copied from original issue: lheagy/simpegem#25

Copied from original issue: simpeg/simpeg#149

utils.spzeros efficiency

I think we want to use sp.dia_matrix instead of sp.csr_matrix for utils.spzeros (2x speedup)

AUTHORS.rst

The author list is out of date. Is this something we want to be maintaining manually? or should we just point to github contributions?

cc @rowanc1

Cylindrical anisotropy

From @rowanc1 on January 5, 2017 19:53

As it is written now, the anisotropy is radial and in z. This is perhaps confusing, especially as theta is necessary but excluded from the tensor.

Copied from original issue: simpeg/simpeg#529

Boundary Conditions

From @lheagy on April 6, 2016 19:18

@sgkang and I were chatting about boundary conditions. Currently we allow two boundary conditions to be set on each boundary, which is likely confusing and potentially dangerous.

Suggestion: in the a general case for a 2D problem something like: (@sgkang feel free to edit this!)

BClist = [('mixed', {'alpha': 0.5, 'beta': 0.5, 'sigma':1}), 'neumann', 'dirichlet', ('neumann',10) ]
Problem(mesh, mapping, BC= BClist) 

Copied from original issue: simpeg/simpeg#286

conda forge

We should get discretize up on conda-forge (https://conda-forge.org/docs/recipe.html). This will help overcome many of the installation issues we keep encountering with windows machines, as described in #32 (cc @sdevriese, @yangdikun).

@leouieda has experience submitting recipes and can hopefully help :)

Crash on install of SimPEG seems to be caused by line in discretize/setup.py

Hi all,

Trying to install SimPEG on my mac. Tried using pip in both python 2.7 and 3.5. I got the same error for both.

Collecting discretize (from SimPEG)
Using cached discretize-0.1.4.tar.gz
Complete output from command python setup.py egg_info:
Traceback (most recent call last):
File "", line 1, in
File "/private/var/folders/kb/40gm7_gx27178vg4v7gswlxh0000gn/T/pip-build-yni1k5cp/discretize/setup.py", line 9, in
from numpy.distutils.core import setup
ImportError: No module named 'numpy'

Looking at the setup.py file I see that on line 9 it says

from numpy.distutils.core import setup

By no means a python expert but looking at other setup.py files and they mainly seem to have the line

from distutils.core import setup

Could this be the problem? If not back to the drawing board. Any help would be appreciated.

Thanks
Sean

mesh.domain

From @lheagy on October 25, 2016 16:44

Often, it is handy to check the extent of a mesh (eg. when making sure it is large enough wrt skin-depths). Right now, I would go through each dimension and print the min and max cell center. Perhaps we can accomplish this through a mesh.domain

@property
domain(self):
    if mesh.dim == 3 # same idea for 1D, 2D 
    return np.array(
        [self.vectorCCx.min(), self.vectorCCy.min(), self.vectorCCz.min()], 
        [self.vectorCCx.max(), self.vectorCCy.max(), self.vectorCCz.max()]
    )

thoughts?

Copied from original issue: simpeg/simpeg#484

Discretize and Properties on Windows Machine

Devi and I are getting the same issue on our Windows machine:

This issue does not seem to arise on Mac.

Properties v0.3.5 (from pip install)
Discretize from dev branch

Thanks!

Test the visualizations

Plotting on meshes

From @lheagy on June 2, 2016 19:48

@thast is visualizing inversion results on SimPEG meshes and has a couple requests

• showPadding = False we should be able to detect where cell sizes change (for tensor at least) and only plot the core region of the mesh (default can be true) - this is important when showing inversion results (there is no data support in the padding). Sometimes you do want to show larger axes to compare lines that are offset but geo-referenced
• mesh.plotCountours same as plotimage - but do contours (or countouf) instead

Copied from original issue: simpeg/simpeg#330

TreeMesh plotslice Colormap

The cmap option needs to be implemented on TreeMesh

PlotSlice

From @rowanc1 on October 27, 2014 20:9

Have the ability to flip the x-y axis:

M.plotSlice(...,flipXY=False,...)

Copied from original issue: simpeg/simpeg#79

OcTree Interpolation

From @rowanc1 on November 25, 2015 17:47

Currently this is O(h), and is nearest neighbour in some of the implementations. This is good for some applications, but likely not for most of the ones that we will actually use interpolation for (Rx, Src, ...).

This needs to be fixed by interpolating over an 6x6 block (4x4 of which will be active). This is a bit involved in the current implementation (and will be more computationally expensive).

Copied from original issue: simpeg/simpeg#182

Mesh IO

From @rowanc1 on January 28, 2016 22:16

I would propose the file structure be replaced:

#----- replace ------
Utils.readUBCTensorMesh('myMesh.msh')
#----- with ------
M = SimPEG.Mesh.TensorMesh.readUBC('myMesh.msh')
v  = M.readModelUBC('model.dat')

This works well with the various types of meshes as well.

M.writeUBC('myNewMesh.msh')
M.writeModelUBC('myModel.dat', v)

I think that would clean up the syntax a bit, and allow these IO formats to live in the correct mesh file, but be used the same.

Similar:
writeVTK

Copied from original issue: simpeg/simpeg#212

Octree model sections

Originally from @micmitch in simpeg/simpeg#575

I have set up a small problem for testing the regularization operators on an octree mesh using a DC problem. This small problem (~50,000 cells) inverts in about 10 min but it takes almost as long for it to generate a 3 panel plot of model sections. I'm assuming this has something to do with interpolation. Any thoughts on how we could potentially speed this up?

pinned sphinx version

in #601, we pinned to sphinx version 1.5.6 as there was a change in the sphinx api that is incompatible with the most recent sphinx gallery. Once sphinx-gallery/sphinx-gallery#241 closes, we can release that and move to the updated sphinx.

Include method gridH to mesh typcs Tensor and Tree

For the code I am writing, I would like to have a method mesh.gridH which behaves like mesh.gridCC. The method would output the cell dimensions of all cells into an array of dimension NX3. I know that from the mesh properties, this can be computed fairly easily. Nevertheless, it would come in handy and avoid cluttering up space in my Problem class.

PEP8

Because I know that we should.

This could probably help ease us in:
http://propertiespy.readthedocs.io/en/latest/content/renamed.html#renamed-property

I think that there are a number of handy properties that people use quite often that would need some thought:

nC, nF, vnE etc. I think true PEP8 doesn't even like the 2 letter variable names? For some of this stuff that is used constantly, I think that is fine?

faceDiv --> face_div

Oh boy.

fatiando/fatiando#278

Thoughts @lheagy?

Fatiando

I am leaving this open for a general open discussion around the idea of bringing some aspects of the fatiando meshing/gridding software to either make use of or replace aspects of this package. With the goal of de-duplicating small parts of the geo-scientific stack and bringing in some common dependencies.

The SquareMesh for example is quite similar to the TensorMesh in 2D:
http://www.fatiando.org/api/mesher.html#fatiando.mesher.SquareMesh

The 3D tensor mesh and octree meshes could generate PrismMeshes:
http://www.fatiando.org/api/mesher.html#fatiando.mesher.PrismMesh

It looks like there would need to be some work in looping over cells to generate geometric primitives. This is actually exactly what we are doing in the treemesh code (this code is dense).
https://github.com/simpeg/discretize/blob/master/discretize/TreeMesh.py#L2282

But basically what I could see is:

mesh = TensorMesh([50, 50])
for square in mesh:
    square.bounds  # fatiando style


mesh = TensorMesh([50, 50, 1])  # What, now in 3D!
for prism in mesh:
    prism.bounds  # fatiando style

There are a couple major differences between the packages, much of the discretize codebase works on vectorization/matrix multiplication to create differential operators. I believe that fatiando is mostly about for loops, and creating (on the fly) representations of the geometric primitives so that calculations can be done on each one?

There may also be some opportunities to work together on how we deal with physical properties, how these are related to the geometric objects, and how you can abstractly define them (e.g. through mapping classes). This was a big push in SimPEG recently, but that got us to v0.9, still lots of thought to be done on how best to do this.

Again, I see this as a way to bring our two communities a bit closer together, as there are significant overlaps: grav/mag, srtomo, inv framework, utilities like data fetching, etc. But one step at a time. :)

@leouieda has noted some of these things in fatiando/fatiando#278 and also pointed out a few opportunities to improve the discretize package (e.g. pep8) #39, which we should seriously consider.

Bug in the 3D cyl mesh interpolation

The size of a face interpolation matrix doesn't match the number of faces

Wheels for Windows

Should consider building wheels for distribution without compilation.

@jcapriot do you have any idea how to do this? I think that we can create a wheel on appveyor, I will give it a first shot.

http://wheel.readthedocs.io/en/latest/
https://packaging.python.org/appveyor/

Extract core mesh from TreeMesh

The plotting for an octree is insanely slow. There is also some times that you want to work with a core mesh as a tensor mesh. This requires you to extract pieces and use those.

A quick sketch of the algorithm in 2D is something like:

# find the first cell of the tree mesh
x0_tensor = MT.gridCC[0, :]
for c in M:
    if np.allclose(c.center, x0_tensor):
        first_cell = c
        break
print(first_cell.center)
print(first_cell._pointer)
ptr = first_cell._pointer

index = np.zeros(MT.vnC, dtype=int)

for ii, nx in enumerate(range(MT.nCx)):
    for jj, ny in enumerate(range(MT.nCy)):
        # if you are not in the last level of the tree mesh 
        # you will have to skip up by that number
        # The last entry of the pointer refers to your level
        ptr_ij = [ptr[0] + nx, ptr[1] + ny, ptr[2]]
        index[ii, jj] = M._cc2i[M._index(ptr_ij)]

index = discretize.utils.mkvc(index)

And then you can plot it:

fig, axs = plt.subplots(1, 2, figsize=(15,8))
v = np.sin(M.gridCC[:, 0] * 2.0 * np.pi) * np.sin(M.gridCC[:, 1] * 2.0 * np.pi) + M.gridCC[:, 1]
M.plotImage(v, ax=axs[0])
axs[0].plot(MT.gridCC[:,0], MT.gridCC[:,1], 'kx')

MT.plotImage(v[index], ax=axs[1])

Would need a bunch of error checking, and extensions to 3D to make this robust.

cc @micmitch

Interpolation for LOM

From @rowanc1 on November 26, 2013 23:18

Need a similar interpolation scheme for the logically orthogonal mesh.

Inverse distance weighting? Delaunay Triangulation?

This will be expensive as the mesh gets big, and in 3D.

---

• Bitbucket: https://bitbucket.org/rcockett/simpeg/issue/23
• Originally Reported By:
• Originally Created At: 2013-11-05T19:41:11.645

Copied from original issue: simpeg/simpeg#25

Expand `vType` in `plot_3d_slicer`

The plot_3d_slicer introduced in v0.3.4 has a vType-argument, but at the moment only CC is implemented.

Expand it to accept same vType-list as plotSlice:
CC, CCv, N, F, E, Fx, Fy, Fz, E, Ex, Ey, Ez.

#110 would resolve this issue. However, I think #110 might be a bit more work, and this issue could be resolve more easily with a few tweaks.

Laplacian operator in DiffOperators

From @rowanc1 on November 26, 2013 23:18

add it

---

• Bitbucket: https://bitbucket.org/rcockett/simpeg/issue/27
• Originally Reported By:
• Originally Created At: 2013-11-06T19:27:56.081

Copied from original issue: simpeg/simpeg#29

Refactoring `plotSlice` for implementation in other plotting routines

A refactoring of plotSlice, in TensorMesh but also in TreeMesh, is required in order to be used by other plotting routines such as plot_3d_slicer.

plotSlice should have an option to only return the actual plot for that axis, without any decorations such as axis labels etc. This could be a pcolormesh instance or a collection of Polygon-patches or another instance. But just the plot, nothing else.

CylMesh Theta Discretization

I think @lheagy has already started on the theta discretization, which is needed for some more accurate wire paths in a casing example.

Later down the road, could be interesting to do both r[0]>0 or \theta \in [0, 1.5]?

Devin Cowan (@dccowan) permission to push branch to discretize

I am contributing towards SimPEG and discretize. I require some permissions.

possible bug in tree mesh readmodelUBC

when reading in a octree model using mesh = Mesh.TreeMesh.readUBC(mesh_file) and then the model using model = mesh.readModelUBC(model_file) it appears that the model is returned as a list of an array, model
Out[62]:
[array([ 3.28385e-03, 3.40776e-03, 5.75453e-03, ..., -9.99900e+03,
-9.99900e+03, -9.99900e+03])]

note extra brackets at start and end.
when reading a tensor model it just returns an array
with no extra brackets

Tab completion on TreeMesh objects

After creating a new TreeMesh object tab completion on other objects/packages works.

Trying to tab complete on the mesh object however, produces no results and breaks any further attempts to tab complete on other objects or packages.

Below is a simple example which shows this.

Here is the mesh file if you would like to try reproducing the problem.

SingleBlk5m_Top_3mFromTunnel_Octree_Core14_10_6_1kmPad_Reg.txt

need simple cellGrad for regularization on TreeMesh

From @micmitch on November 11, 2016 3:25

TreeMesh does not inherit the functions in DiffOperators.py this currently limits us to the use of BaseRegularization which does not use any derivative terms for model smoothness. Hopefully we can use the faceDiv operator that we have to create a simple gradient operator for use in the regularization. Just not sure where this should be built or where to start from.

Copied from original issue: simpeg/simpeg#504

cache the averaging matrices

right now, we use sp.hstack on every call to the averaging matrices to stack the contributions from the x, y, z components (if in 3D). This can be a slow operation for large meshes, so we should use lazy loading to cache the averaging matrices.

See simpeg/simpeg#587

plotting without padding cells

mesh.plotImage(sigma, padding=False)

would be very useful and help speed up a lot of the plotting code

Import DC Octree Mesh from UBC fails

From @thast on July 27, 2016 0:22

@rowanc1 : is it your function? or @sgkang ?

Copied from original issue: simpeg/simpeg#368

Bug in creating and/or plotting Tree meshes

From @dccowan on August 15, 2017 20:53

I have noticed that the x0 position in Tree mesh objects is not respected. As a result, some of the following occurs:

• You can create a tree mesh object where x0 isn't the origin. However when you want to plotGrid(), it ALWAYS plots with x0=(0,0,0). This is strange given that you can use meshObj.x0 to verify that the bottom southwest corner is not at the origin.
• If you want to refine the mesh using meshObj.refine(), it doesn't take into account the x0 position.
• When you want to create a model, it DOES take into account the x0 position. Which means the model looks right, but gets shifted over when you plot it.

I have a feeling that this is a pretty easy fix. I just don't know where to find the specific lines of code.

Copied from original issue: simpeg/simpeg#632

octree mesh reshape (mesh.r)

Given the recent work on octree there are a couple of tensor mesh functions that would be useful to implement in octree. One is mesh.r . From the tensor mesh

r is a quick reshape command that will do the best it can at giving you what you want.
For example, you have a face variable, and you want the x component of it reshaped to a 3D matrix.
r can fulfil your dreams:

I would suggest using a better name such as reshape?

One way to do it might be to resample back to Tensor using something like this (thanks Dom)

Can my dreams be fulfilled?
`
if isinstance(meshInput, Mesh.TensorMesh):
tree = cKDTree(mesh.gridCC)

    dd, ind = tree.query(meshInput.gridCC) # can load any Tensor mesh here.

    # Model out

    vec_xyz = Utils.matutils.atp2xyz(mrec_MVI_S.reshape((nC, 3), order='F')).reshape((nC, 3), order='F')
    vec_x = actvMap * vec_xyz[:,0]
    vec_y = actvMap * vec_xyz[:,1]
    vec_z = actvMap * vec_xyz[:,2]

    vec_xyzTensor = np.zeros((meshInput.nC,3))
    vec_xyzTensor = np.c_[vec_x[ind], vec_y[ind], vec_z[ind]]

    PF.MagneticsDriver.writeVectorUBC(meshInput, work_dir+out_dir + 'MVI_S_Tensor.fld', vec_xyzTensor)

    amp = np.sum(vec_xyzTensor**2., axis=1)**0.5
    meshInput.writeModelUBC(work_dir+out_dir + 'MVI_S_Tensor.amp', amp)
    meshInput.writeUBC(work_dir+out_dir + driver.mshfile.split('/')[2])

`

Installing from pypi

Currently, from pypi, the build_ext --inplace is not being run, so efficiency warnings are thrown.

For example, when testing SimPEG https://travis-ci.org/simpeg/simpeg/jobs/201040464

Bug in travis early exits

The latest pypi release deployed earlier than anticipated (thanks @jcapriot for catching this). It wasn't a problem as the deployed version was what we intended to deploy - however, it should have been deployed only upon the merge to master or the tagged release...

after_success:

  - bash <(curl -s https://codecov.io/bash)

  # early exit if not on a deplotment branch
  - if ! [ "$TRAVIS_BRANCH" = "$MASTER_BRANCH" -o "$TRAVIS_TAG" = "true" ]; then
      echo "Not deploying (because this is not a deployment branch)" ;
      exit 0 ;
    fi
  - if [ "$TRAVIS_PULL_REQUEST" = "true" ]; then
      echo "Not deploying (because this is a pull request)" ;
      exit 0 ;
    fi

pull requests directly onto dev

Similar to simpeg/simpeg#651: To speed up getting features live into discretize, we could do pull requests directly onto master.

Impacts

• your features are rapidly live to the world!
• pull requests should be small and focussed
• pull request notes will be used as the release notes and so should contain sufficient detail to describe the changes, any breaks in backwards compatibility, etc.

TODO's

• switch main branch to master on github
• ensure dev and master synced
• change base of current pr's onto dev to master
• delete dev branch
• update travis to remove deployment sequence from dev

Timeline

We can do this as soon as #73 is in. Are there any concerns about doing this soon? cc @rowanc1, @grosenkj, @sgkang, @thast, @fourndo. Can you please give a 👍 or 👎?

TensorMesh object has no attribute 'aveCC2FV'

TensorMesh used to have this attribute (aveCC2FV), late 2017, but the latest version seems to have dropped it.
Lindsey suggested :

In the mean time, I believe ‘aveF2CCV’ is still there. You could use ‘aveCCV2F = Utica.sdiag(1/mesh.aveF2CCV.T.sum(1)) * mesh.aveF2CCV.T’

The sdiag term will make sure that the ends of the mesh are dealt with (eg not divided by 2)

but this is pretty cumbersome when a previous attribute already did it.

Could aveCC2FV please be added back? I'm using it in an SP inverison

Mesh Tensor Generation

From @rowanc1 on April 6, 2016 0:29

@lheagy was saying that the mesh tensor needs a bit more work, would be good to do something like:

[ ( 0.05, {"to":100}), (0.05, {"max":25}, 1.3), (25, {"to":1000}) ]

Thoughts on the API?

Copied from original issue: simpeg/simpeg#283

Separating operator stencils and geometry

Right now, there is a mix of the existence of Stencil operators, which contain the core of the operator, but not the length scales associated with a specific mesh, and those that don't.

The stencil is specific to the mesh type, but the way we incorporate sizes (and do deflation) is general and could be abstracted. This would clarify what pieces need to be implemented when designing a new mesh and which can be inherited. Also, components can (and should!) be broken up so that when calculating the Div for example, it just grabs the x, y, z div pieces (this means they can then be tested and used independently.

This would require a bit of re-organization, but should make it more clear which code needs to live with each specific mesh, and what comes from DiffOperators

cc. @rowanc1

Tree Mesh on Windows

This currently does not work. If someone is interested in helping, it has something to do with long integers that are not represented well in c?

The cython is doing some bitshifting which windows does not like... :(

docs are out of date

The discretize docs last built for version (e.g.: http://discretize.simpeg.xyz/en/latest/auto_examples/plot_quadtree_divergence.html)

getEdgeInnerProduct invMat = True needs to be implemented

originally from simpeg/simpeg#105

Should be able to call: (for diagonal models, it is simple, for fully anisotropic, it involves a matrix solve)
self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u)
self.mesh.getFaceInnerProductDeriv(self.curModel.mui, invMat=True)(u)

Interpolation

From @rowanc1 on November 17, 2015 22:0

In the interp utils we should be using np.searchsorted rather than np.argmin over the vectors:

Speed up of O(log(n)) vs O(n) for each point you are interpolating, could be significant, especially on large meshes.

Copied from original issue: simpeg/simpeg#170

inner product names

• getEdgeInnerProduct --> edgeInnerProduct
• getFaceInnerProduct --> faceInnerProduct

thoughts @rowanc1 ?

meshutils.exampleLrmGrid --> examples ?

originally from: simpeg/simpeg#278

this util code looks better suited for Examples to me, thoughts @rowanc1 ?

Concave Areas not calculated properly in LOM

From @rowanc1 on November 26, 2013 23:17

area calculation (i.e. 'volume' in 2D) is not correct if the cell is concave.

Currently averaging the four parallelograms, and I think this works out to give the area of the convex hull of the cell.

Need to write test cases and explore this issue.

---

• Bitbucket: https://bitbucket.org/rcockett/simpeg/issue/3
• Originally Reported By:
• Originally Created At: 2013-08-01T17:33:02.126

Copied from original issue: simpeg/simpeg#5