The branch extended_testing contains a script test/regression/check_tutorials.py that runs all tutorials (e.g. all *.ipynb in the tutorials folder) and leaves an assertion error if something goes wrong. And something indeed goes wrong.

@alessiofumagalli @rbe051 @IvarStefansson Please check that 'your' tutorials are working, or make them work. When this is done, and I've fixed mine I will try to add this job to the daily run.

Minimal example:

`
import numpy as np

from porepy.grids import structured
from porepy.params import fourth_order_tensor, bc
from porepy.numerics.fv import mpsa

n = 5
g = structured.CartGrid([n,n])
g.compute_geometry()

dirich = np.ravel(np.argwhere(g.face_centers[1] < 1e-10))
bound = bc.BoundaryCondition(g, dirich, ['dir']*dirich.size)

lam = np.ones(g.num_cells)
mu = np.ones(g.num_cells)
constit = fourth_order_tensor.FourthOrderTensor(g.dim, mu, lam)

_, _ = mpsa.mpsa(g, constit, bound)

_, _ = mpsa.mpsa(g, constit, bound)
`

Should be relatively straightforward, but should be done.

While fractures that intersect faces of the bounding box can be truncated, it is tacitly assumed that no fractures intersects the edges of the box. This needs to be fixed, it is e.g. relevant when compartments in the domain (say, aquifers, aquitards etc) are imposed as fractures.

A fix will require a rewrite of how the domain is exported to gmsh (in the future also other software packages).

indPtr = np.array([1,5,5,6])
select = np.array([0,1,2])
mcolon(indPtr[select], indPtr[select + 1]) # OK

select = np.array([0,1])
mcolon(indPtr[select], indPtr[select + 1]) # Crashes

Traceback (most recent call last):
File "", line 1, in
File "/home/runar/uib/porepy/src/porepy/utils/mcolon.py", line 59, in mcolon
print(n, 'n')

It seems like the intersection algorithm fails when a fracture vertex lie on a fracture edge.
It is something with the ordering of the nodes. If the fracture has only three vertices it is ok.

import numpy as np
from porepy.fracs.fractures import Fracture
from porepy.fracs import meshing
f_1 = Fracture(np.array([[0,0,0], [1,0,0],[1,1,0], [0,1,0]]).T)
f_2 = Fracture(np.array([[0,0,1], [0,0,-1], [1,1,-1]]).T)
gb = meshing.simplex_grid([f_1, f_2])
'

Identification of fractures is now done by the fracture index, which is assigned in an as systematic way as possible, but is still not robust. Should be fixed.

Mappings from higher and lower dimensional grids both to edges.

Replace ad hoc print statements with proper logging, using standard Python tools.

Uniform the how vem handle the scalar source term with finite volume, where the source term is considered already integrated.

we now have:
mcolon(np.array([0]), np.array([2]))
[0,1,2]

however, it would be more "pythonic" if it returns
[0,1]

as np.arange does:
np.arange(0,2)
[0,1]

In mpsa._zero_neu_rows, line 1253, I run out of memory for a problem that ought not cause problems. Likely reason is that the combination of sparse matrices do not work as expected.

Minimal example (32GB RAM):
n=30
g = structured.CartGrid([n, n, n])
g.compute_geometry()
k = fourth_order_tensor.FourthOrderTensor(g.dim, mu=np.ones(g.num_cells),
lmbda=np.ones(g.num_cells))
bound_faces = g.get_boundary_faces()
dir_bound = np.argwhere(g.face_centers[1, :] < 1e-10)
bound = bc.BoundaryCondition(g, dir_bound.ravel(), ['dir'] * dir_bound.size)
stress, bound_stress = mpsa.mpsa(g, k, bound, inverter='numba')

Implement a data class as the main (ideally only) means to pass around parameters etc. Replaces #24 , extends #25 .

test.unit.test_dual_vem.BasicsTest.test_dual_vem_2d_ani_cart_surf() fails occasionally, with no clear reasons. See e.g. Travis build 482.

Apparently, the test alternates between two different answers, with one being the identified correct one.

Should not this work?

f = np.array([[0,0,0.5], [0,1,0.5], [1,1,0.5], [1,0,0.5]]).T
gb = meshing.simplex_grid([f])

1. It should be possible to save a mesh.
2. Save simulation result, in a way that enables restart.
3. Save Parameter class, or its extension.

@IvarStefansson @rbe051 The emerging structure, where the discretization is provided by module-level functions, but also with a class wrapper, is confusing. I think we do want to keep the classes, both to be compatible with the approach in VEM, and also for convenient handling of mixed-dimensional problems. Also note that the classes are shells that do not store data (thus they are quite equivalent with functions, but appears more elegant to me).

Do you see any downside with dropping the functions, and leave a single, class, entry-point to the discretizations?

I just did a clean install of porepy and conda, but got a issue when trying to install the requirements in pip.

I did:
conda create --name porepy python=3.6
source activate porepy
pip install -r requirements-dev.txt

However, something broke in the downloading of pymetis:
ModuleNotFoundError: No module named 'six'

I installed six with pip, and after that everything ran smoothly

Needed for mortars.

Is it possible to import two m x n arrays depicting permeability and porosity, respectively, and use that to solve for single-phase flow in porepy? I would like to solve for flux (flow rate) with given boundary conditions on a Cartesian mesh that is simply the imported array without any discrete fracture because the imported permeability (and porosity) values are implicitly accounting for the presence of fractures. Another related question - can I use the pressure at each time step to modify my permeability array and use that modified permeability for the pressure/flux calculation on next time step? It would be really useful to have this capability in porepy if it's not there currently.

I'll submit an update of the exporter soon.
The idea is to introduce a class structure and for grids which are not changing in time store them (compressed as suggested in https://stackoverflow.com/questions/19500530/compress-python-object-in-memory) as vtk grid object. This is an option and time dependent grid will be handled as well. For big mesh the compromise between create the vtk grid each time step or store it need a future investigation.

I was trying to run an example, but it seems like none of the given examples is working for me because of ModuleNotFoundError. Specifically, import porepy works for me (Python 3), but any other import from porepy as used in examples do not work for me. I have tested my installation (from source) and it passed the testing. For instance, none of the following imports used in examples work for my installation:

from porepy.viz import exporter
from porepy.fracs import importer
from porepy.params import bc, tensor
from porepy.params.data import Parameters
from porepy.numerics.mixed_dim import coupler
from porepy.numerics.fv import tpfa, mpfa
from porepy.utils.errors import error

Sorry I cannot provide a simple example, see the branch update_coarsening.

It seams that with different version of python the map between face (of higher dimension) to cell (of lower dimension) changes during the runs.

Moreover with version 3.5 and 2.7 the result is not stable, sometimes the tests work and sometimes not. So far, only 3.6 is stable with that respect.

This is a major problem with testing for different python versions.

The L intersections are not working in 2d. Minimal example:

from porepy.viz import plot_grid
from porepy.fracs import importer

mesh_kwargs = {}
mesh_kwargs['mesh_size'] = {'mode': 'constant',
                            'value': 0.5, 'bound_value': 0.5}

domain = {'xmin': 0, 'xmax': 1, 'ymin': 0, 'ymax': 1}
gb = importer.from_csv('single.csv', mesh_kwargs, domain)
gb.compute_geometry()
plot_grid.plot_grid(gb, info="all", alpha=0)

where single.csv is

FID,START_X,START_Y,END_X,END_Y
1,0.2,0.5,0.8,0.5
2,0.2,0.1,0.2,0.5

I get the following error

    gb = importer.from_csv('single.csv', mesh_kwargs, domain)
  File "/home/elle/ANIGMA/porepy/src/porepy/fracs/importer.py", line 56, in from_csv
    return meshing.simplex_grid(f_set, domain, **mesh_kwargs)
  File "/home/elle/ANIGMA/porepy/src/porepy/fracs/meshing.py", line 72, in simplex_grid
    grids = simplex.triangle_grid(frac_dic, domain, **kwargs)
  File "/home/elle/ANIGMA/porepy/src/porepy/fracs/simplex.py", line 253, in triangle_grid
    pts, cells, phys_names, cell_info, line_tag=const.PHYSICAL_NAME_FRACTURES)
  File "/home/elle/ANIGMA/porepy/src/porepy/grids/gmsh/mesh_2_grid.py", line 157, in create_1d_grids
    g = create_embedded_line_grid(loc_coord, loc_pts_1d)
  File "/home/elle/ANIGMA/porepy/src/porepy/grids/gmsh/mesh_2_grid.py", line 186, in create_embedded_line_grid
    tangent = cg.compute_tangent(loc_coord)
  File "/home/elle/ANIGMA/porepy/src/porepy/utils/comp_geom.py", line 1560, in compute_tangent
    assert not np.allclose(tangent, np.zeros(3))
AssertionError

Probably has to do with recent updates of physical names.

Until fixed, meshio 1.9 should be safe.

I was trying to clone the project on Windows and got this invalid argument error.
@alessiofumagalli can you please change names to folders and replace ':' perhaps with '_'?

GridBucket has suffered from an inflation of methods, with very similar names to boot. Case in point:
node_prop(), node_props(), nodes_prop(), similar for edge_*.

I will clean this up, at least by moving to more descriptive and distinct names, and possibly also by killing some methods. @alessiofumagalli @IvarStefansson @rbe051 : Can you please guide me to where these methods are used, so that I can get an idea of their usage?

We should also improve test coverage for the bucket. I will start on this, but will likely need inspiration from others.

I came over a case where a face in the 3D grid is missing. The face should be a fracture face. The face exist in the 2D grid (as a cell), however, the 2D cell does not have a mapping to any face in the 3D domain. I have checked that it is not any problem with the mapping, and indeed the 3D face does not exists. I assume this also means that the two cells that should be connected to the face (note that this is before we split the grid) are also missing.

Either, we are not reading the output from gmsh improperly, or gmsh removes this cell for some reason...

import numpy as np

from porepy.fracs import fractures, meshing
from porepy.viz.exporter import export_vtk


def create_fractures(data):
    data = np.atleast_2d(data)
    centers = data[:, 0:3]
    maj_ax = data[:, 3]
    min_ax = data[:, 4]
    maj_ax_ang = data[:, 5]
    strike_ang = data[:, 6]
    dip_ang = data[:, 7]
    if data.shape[1] == 9:
        num_points = data[:, 8]
    else:
        num_points = 16 * np.ones(data.shape[0])

    frac_list = []

    for i in range(maj_ax.shape[0]):
        frac_list.append(fractures.EllipticFracture(centers[i, :],
                                                    maj_ax[i],
                                                    min_ax[i],
                                                    maj_ax_ang[i],
                                                    strike_ang[i],
                                                    dip_ang[i],
                                                    num_points[i]))
    fracture_network = fractures.FractureNetwork(frac_list)
    return fracture_network


frac_list = np.array([[2000, 2000, 2000, 1500, 1500, 0, 100, 80],
                      [1800, 1800, 2150, 700, 700, 0, 170, 70],
                      [2650, 2200, 1700, 850, 850, 0, 130, 70],
                      [2100, 2100, 2000, 1000, 1000, 0, 60, 80]])

box = {'xmin': 0, 'ymin': 0, 'zmin': 0,
       'xmax': 4000, 'ymax':  4000, 'zmax': 5000}
f = create_fractures(frac_list)
g = meshing.simplex_grid(f, box)
g.assign_node_ordering()

export_vtk(g, 'temp')
`

In fractured domains, in the construction of the individual grids, we have accessible the mapping between 2d grids and fracture surfaces as specified by users. However we do not make use of the information. As a result, when individual grids / cells on fracture surfaces are accessed later, e.g. for data assignment, cumbersome and error prone routines must be followed (often based on geometry).

Proposed solution: It should be possible to inquire the grid_bucket for getting grids based on the fracture id.

A similar problem exist for 2d domains as well.

The cart_grid_2d function fails if there is more than 1 intersection. Minimal example:

    frac1 = np.array([[1, 4], [2, 2]])
    frac2 = np.array([[2, 2], [1, 4]])
    frac3 = np.array([[3, 3], [1, 4]])
    fracs = [frac1, frac2, frac3]
    gb = porepy.fracs.structured.cart_grid_2d(fracs, [5, 5])

Solution: Line 229,

for global_node in np.where(shared_nodes > 1):

Should be changed to

for global_node in np.where(shared_nodes > 1)[0]:

(I could have fixed it myself if I had access...)

What do you think? We basically use only embedded grids now and keeping track of it can be annoying. Also the True option is always correct (maybe takes a bit more time) the other no.

Boundaries of parameter subdomains (e.g. aquifer/aquitard, under and overburden) are conveniently introduced by fracture planes. However, these can currently not extend to the edges of the bounding box, as the latter is tacitly assumed not to be intersected by other objects.

The modification will require reconfiguration of the way the boundary is introduced in the meshing.

For large problems, the split_fractures part of the meshing algorithm takes exceedingly long time. Below follows run script an output to illustrate the problems.

Run with branch frac_mesh_3d, commit e765a57 .

Run script

import numpy as np

from porepy.fracs import extrusion, importer, meshing
from porepy.fracs.fractures import Fracture, FractureNetwork
from porepy.viz import exporter

np.random.seed(7)
pt, edges = importer.fractures_from_csv('barely-connected_4pp.log
barely-connected_4pp.log
')
fractures = extrusion.fractures_from_outcrop(pt, edges)

f = fractures

network = FractureNetwork(f)
network.tol = 1e-5
network.find_intersections()
network.split_intersections()

gb = meshing.simplex_grid(network=network, h_ideal=0.5, h_min=0.1, verbose=1, ensure_matching_face_cell=False)

########## Output ########################

Construct mesh
Use existing intersections
Use existing decomposition
Gmsh processed file successfully

Grid creation completed. Elapsed time 1.3945231437683105

Created 1 3-d grids with 81662 cells
Created 66 2-d grids with 16949 cells
Created 82 1-d grids with 234 cells

Done. Elapsed time 75.7353572845459
Assemble in bucket
/home/eke001/Dropbox/workspace/python/porepy/src/porepy/fracs/meshing.py:404: UserWarning: Found inconsistency between cells and higher
dimensional faces. Continuing, faces crossed
dimensional faces. Continuing, faces crossed''')
Done. Elapsed time 10.409296751022339
Compute geometry
Done. Elapsed time 0.6946992874145508
Split fractures
Done. Elapsed time 497.07970118522644
Mesh construction completed. Total time 583.9497916698456

Hi, I have a problem when I couple the extrusion and the non-conforming algorithm. However, I think the main problem is due to the latter. Note: if considering the extrusion couple with a conforming algorithm it works.
Minimal example:

import numpy as np                                                               
import scipy.sparse as sps                                                       
                                                                                 
from porepy.viz.exporter import Exporter                                         
from porepy.fracs import importer                                                
                                                                                 
from porepy.fracs import extrusion                                               
from porepy.fracs.fractures import Fracture, FractureNetwork                     
from porepy.fracs import meshing                                                 
                                                                                                                                                               
folder_export = 'example_4_vem/'                                                 
file_export = 'vem'                                                              
tol = 1e-8                                                                       
                                                                                 
# The fractures are not aligned with the mesh                                    
pts = np.array([[0, 0], [2, 2], [1, 1], [1, 3]]).T                               
edges = np.array([[0, 1], [2, 3]]).T                                             
                                                                                 
extrustion_kwargs = {}                                                           
extrustion_kwargs['tol'] = tol                                                   
extrustion_kwargs['exposure_angle'] = 2*np.pi/3.*np.ones(edges.shape[1])         
                                                                                 
fractures = extrusion.fractures_from_outcrop(pts, edges, **extrustion_kwargs)    
network = FractureNetwork(fractures, tol=tol)                                    
                                                                                 
# If uncomment the code work with the new network                                
# f_1 = Fracture(np.array([[-1, -1, 0], [1, -1, 0], [1, 1, 0], [-1, 1, 0]]).T)   
# f_2 = Fracture(np.array([[-1, 0, -1], [1, 0, -1], [1, 0, 1], [-1, 0, 1]]).T)   
# f_3 = Fracture(np.array([[-0.75, 0.25, -0.75], [0.75, 0.25, -0.75], [0.75,     
#                            0.25, 0.75], [-0.75, 0.25, 0.75]]).T)               
# network = FractureNetwork([f_1, f_2, f_3])                                     
                                                                                 
mesh_kwargs = {}                                                                 
mesh_kwargs['mesh_size'] = {'mode': 'weighted',                                  
                            'value': 0.5,                                        
                            'bound_value': 0.5,                                  
                            'tol': tol}                                          
                                                                                 
gb = meshing.dfn(network, conforming=False, **mesh_kwargs)                       
                                                                                 
gb.remove_nodes(lambda g: g.dim == 0)                                            
gb.compute_geometry()                                                            
gb.assign_node_ordering()                                                        
                                                                                 
save = Exporter(gb, file_export, folder_export)                                  
save.write_vtk()

Special case of T intersection of 2d planes. The intersection line is interpreted as a "boundary line" of the base of the T, possibly because it reaches the boundary of the other plane. This ultimately yields
assert sorted_lines[0, 0] == sorted_lines[1, -1]
AssertionError
Happens for f_1 with both f_2 and f_3. but not f_4 (corner to corner):

def test_issue_56():

domain = {'xmin': -2, 'xmax': 2, 'ymin': -2, 'ymax': 2, 
          'zmin': -2, 'zmax': 2} 
mesh_size = {'mode': 'constant', 'value': .4, 'bound_value': 1} 
kwargs = {'mesh_size': mesh_size, 'return_expected': True} 
 
f_1 = np.array([[0, 1, 1, 0], [0, 0, 1, 1], [0, 0, 0, 0]]) 
f_2 = np.array([[0, 1, 1, 0], [0, .5, .5, 0], [0, 0, 1, 1]]) 
f_3 = np.array([[0, 1, 1, 0], [0.3, .5, .5, 0.3], [0, 0, 1, 1]]) 
f_4 = np.array([[0, 1, 1, 0], [0, 1, 1, 0], [0, 0, 1, 1]]) 
grids = meshing.simplex_grid([f_1, f_2], domain,**kwargs) 

In the case of single grid, at the end of the constructor this function is called

update_boundary_face_tag()

which updates only the boundary faces with the tag BOUNDARY. However, in the case of grid bucket the faces at the boundary are also DOMAIN_BOUNDARY, while the BOUNDARY includes also the FRACTURE faces. Can I simply change the line in the previous function from:

self.add_face_tag(bd_faces, FaceTag.BOUNDARY )

to

self.add_face_tag(bd_faces, FaceTag.BOUNDARY | FaceTag.DOMAIN_BOUNDARY)

I installed porepy (from source using Python 2.7) and then followed-up by testing the build as suggested using the following commands:

pip install -r requirements-dev.txt
python setup.py test

However, the test gives the following error:

==================================== ERRORS ====================================
___ ERROR collecting examples/example8/test_advection_diffusion_coupling.py ____
ImportError while importing test module '/home/gary/porepy/examples/example8/test_advection_diffusion_coupling.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
examples/example8/test_advection_diffusion_coupling.py:139: in <module>
    gb = importer.from_csv(folder + 'network.csv', mesh_kwargs, domain)
src/porepy/fracs/importer.py:43: in from_csv
    return meshing.simplex_grid(f_set, domain, **mesh_kwargs)
src/porepy/fracs/meshing.py:90: in simplex_grid
    grids = simplex.triangle_grid(frac_dic, domain, **kwargs)
src/porepy/fracs/simplex.py:270: in triangle_grid
    return triangle_grid_from_gmsh(file_name, **kwargs)
src/porepy/fracs/simplex.py:288: in triangle_grid_from_gmsh
    **gmsh_opts)
src/porepy/grids/gmsh/gmsh_interface.py:469: in run_gmsh
    config = read_config.read()
src/porepy/utils/read_config.py:46: in read
    raise ImportError(s)
E   ImportError: Could not load configuration file for PorePy.
E   To see instructions on how to generate this file, confer help
E   for this module (for instance, type 
E     porepy.utils.read_config? 
E   in ipython)
------------------------------- Captured stdout --------------------------------
/home/gary/porepy/examples/example8/
============================ pytest-warning summary ============================
WC2 /home/gary/porepy/test/unit/test_doctests.py cannot collect 'test_suite' because it is not a function.
!!!!!!!!!!!!!!!!!!! Interrupted: 1 errors during collection !!!!!!!!!!!!!!!!!!!!
================= 1 pytest-warnings, 1 error in 139.03 seconds =================

minimal example:

import numpy as np

from porepy.fracs import fractures, meshing
from porepy.utils import comp_geom

num_fracs = 1
np.random.seed(1)
centers = np.random.rand(3, num_fracs) + 1.5
maj_ax = 0.5 + np.random.rand(num_fracs)
min_ax = maj_ax
rot = np.pi * np.random.rand(num_fracs)
strike = np.pi * np.random.rand(num_fracs)
dip = np.pi * np.random.rand(num_fracs)
frac_list = []
for i in range(num_fracs):
    frac_list.append(fractures.EllipticFracture(centers[:, i],
                                                maj_ax[i],
                                                min_ax[i],
                                                rot[i],
                                                strike[i],
                                                dip[i]))

box = {'xmin': 0, 'ymin': 0, 'zmin': 0,
       'xmax': 4, 'ymax': 4, 'zmax': 4}
gb = meshing.simplex_grid(frac_list, box)
gb.assign_node_ordering()
gb.compute_geometry()

gb.add_node_props(['param'])
for g, d in gb:
    if g.dim == 2:
        comp_geom.map_grid(g)

For certain fracture networks, an error occurs in the function duplicate_nodes(), around L480. The symptom
is indexing with an index iv that is too long compare to the array length.

It seems the bug was introduced in d34fff4, or one of the preceding updates.

I have not managed to recreate the error on a small network (<40 fractures 3D), let me know when you need a case.

I'm preparing a jupyter notebook to solve a darcy and transport in a dfm setting. Trying to be as user friendly as possible.

Viscosity should be added to the parameter class. Also need the figure out the best way to include in the equations.

Below included example on frac_mesh_3d leads to intersecting face pairs. Gmsh says

!! Found 108 pairs of faces are intersecting.

Self-intersection seconds: 0.113017
Writing nodes.
Writing faces.
Error : ------------------------------
Error : Mesh generation error summary
Error : 0 warnings
Error : 1 error
Error : Check the full log for details
Error : ------------------------------
Gmsh failed with status 256

CODE

from porepy.fracs import meshing
import numpy as np
from porepy.viz.exporter import export_vtk

def define_grid():
"""
One horizontal and one vertical 1d fracture in a 2d matrix domain.
"""
mesh_kwargs = {}
mesh_size = 5e-1
mesh_kwargs['mesh_size'] = {'mode': 'constant',
'value': mesh_size, 'bound_value': 1.2*mesh_size}
mesh_kwargs['file_name'] = 'bounding_box_test'

domain = {'xmin': 0, 'xmax': 1,
          'ymin': 0, 'ymax': 1,
          'zmin': 0, 'zmax': 1}
f_1 = np.array([[.5,.5,.5,.5],
                [.4,.5,.5,.4],
                [0.2,0.2,.8,.8]])
f_2 = np.array([[0,.8,.8,0],
                [.5,.5,.5,.5],
                [0.2,0.2,.8,.8]])#[0,0,1,1]])


fracs = [f_1,f_2]
gb = meshing.simplex_grid( fracs,domain,**mesh_kwargs)
gb.assign_node_ordering()

return gb

gb = define_grid()
export_vtk(gb, 'bounding_box_test_grid')

from porepy.viz import exporter                                       
from porepy.fracs import importer

mesh_kwargs = {}                                                                 
mesh_kwargs['mesh_size'] = {'mode': 'constant', 'value': 0.25, 'bound_value': 2} 
gb, _ = importer.from_csv("single.csv", mesh_kwargs)
exporter.export_vtk(gb, 'grid')

where the csv file is

FID,START_X,START_Y,END_X,END_Y                                                  
1,0.5,0,0.5,1                                                                    
3,0.75,0.5,0.75,1                                                                
4,0.49,0.625,0.75,0.625                                                          
5,0.49,0.62,0.75,0.62

Should be applicable to 'all' conservative numerical methods for flow. Formulate as a general DG method, but only implement lowest order for now.

I was wondering if you can help guide how to change the fluid flow physics in the code. I am not an expert in numerical methods and discretizations, so I am wondering if there's an easy way to add/modify a partial differential equation describing some physical phenomenon? Thanks.

Networkx seems to have changed the underlying data structure when updating to 2.0. We're fixed at <2 for now, but need to move on at some point. Or write our own graph class.

`frac1 = np.array([[1,2,4], [1,4,1], [2,2,2]])

frac2 = np.array([[2,2,2], [2,4,1], [1,2,4]])

fracs = [frac1, frac2]

domain = {'xmin': 0, 'ymin': 0, 'zmin': 0,
'xmax': 5, 'ymax': 5, 'zmax': 5,}

path_to_gmsh = .... # Set the sytem path to gmsh

gb = tetrahedral_grid(fracs, domain, gmsh_path = path_to_gms)`

which gives following error:
Find intersections between fractures
Traceback (most recent call last):
File "", line 1, in
File "/home/runar/uib/porepy/src/porepy/fracs/simplex.py", line 92, in tetrahedral_grid
if not network.has_checked_intersections:
File "/home/runar/uib/porepy/src/porepy/fracs/fractures.py", line 869, in find_intersections
isect, _, _ = first.intersects(second, self.tol)
File "/home/runar/uib/porepy/src/porepy/fracs/fractures.py", line 345, in intersects
assert u_bound_pt.shape[1] == 1
IndexError: tuple index out of range

Motivated by comment from @rbe051, I suggest the following renaming:
from_csv() -> mesh_from_csv()
fractures_from_csv() -> lines_from_csv().

Opinions? If not, I will do the restructuring.

Mostly needed for visualizations for FV methods.