Negative number of tetrahedra in high-res mesh generation

Question

Hello,

I am trying to generate a high resolution mesh (dx = 2 um) with mesher, but the terminal output shows a negative number of tetrahedra. When I run it in a simulation, I get "std::bad_array_new_length" errors. Is there anything that can be done to get around this? 

This is my mesher script:

/u1/jav630/Jan10_openCARP-Actors/openCARP-Actors/environment_setup/openCARP/_build/bin/mesher \

-size[0] 0.1444 \

-size[1] 0.2604 \

-size[2] 0.159 \

-bath[0] -0.0 \

-bath[1] -0.0 \

-bath[2] -0.0 \

-center[0] 0.0722 \

-center[1] 0.1302 \

-center[2] 0.0795 \

-resolution[0] 2.0 \

-resolution[1] 2.0 \

-resolution[2] 2.0 \

-mesh 2024-08-01_Fine/block \

-Elem3D 0 \

-fibers.rotEndo 0.0 \

-fibers.rotEpi 0.0 \

-fibers.sheetEndo 90.0 \

-fibers.sheetEpi 90.0 \

-periodic 0 \

-periodic_tag 1234 \

-perturb 0.0

Please let me know if further information is required.

Thank you!

Joyce 

Answer 1

Mesher is a rather simple tool for generating small meshes on the fly. I don't think it's suitable for generating meshes with >3 billion elements.

https://gmsh.info might be a good option here.

Comments

Thanks for getting back to me. I have been trying to use Gmsh, but it doesn't seem to integrate well with openCARP for me. When I convert the .msh file to the openCARP format using meshtool, then run it with openCARP, I keep getting "L5 : error: Empty stimulus[0] electrode def! Aborting! ", even though the same code works with other (Mesher) meshes. Are there any examples you can point me to of how to convert between Gmsh and openCARP meshes? Please let me know if you'd like to see any of the files I am using.

---

How do you define the stimuli? The carputils convenience functions
block_bc_opencarp() / block_boundary_condition() only work with on-the-fly generated Block meshes:
https://opencarp.org/doxygen-carputils/master/namespacecarputils_1_1mesh_1_1block.html

See openCARPs Stimulation -> Stim -> elec documentation for details how to set up stimuli on arbitrary meshes

---

Previously I was defining the stimulus as something like

num_stim = 1
stimulus[0].name = "S1"
stimulus[0].stimtype = 0
stimulus[0].strength = 0.0
stimulus[0].duration = 2.0
stimulus[0].start    = 0.
stimulus[0].npls     = 1
stimulus[0].x0 = 0.0
stimulus[0].xd = 222.0
stimulus[0].y0 = 0.0
stimulus[0].yd = 222.0
stimulus[0].z0 = 0.0
stimulus[0].zd = 222.0

Then in response to your comment I changed it to something like
num_stim = 1
stim[0].name = "S1"
stim[0].elec.geom_type = 2
stim[0].elec.p1 = 222.0
stim[0].elec.p0 = 0.0
stim[0].ptcl.duration = 2.0
stim[0].pulse.strength = 300

but unfortunately I am still getting the same error ("L5 : error: Empty stimulus[0] electrode def! Aborting!"), even though I have not defined any "stimulus" objects. Do you know what I am doing wrong here?

Also, I'm not sure if this makes a difference, but I am using a parameter file (s1.par) and running it with "openCARP +F s1.par", rather than using a Python script.

---

Please provide a minimum example to reproduce this issue

---

I am not sure the best way to send it as I can't attach anything to comments, but I will try and paste the text here. Please let me know if anything is unclear.

This is my Gmsh file, bidomain_005cm3.geo , including instructions for how to run everything at the top :

// -----------------------------------------------------------------------------
//
//  This Gmsh script attempts to produce a simple cuboidal mesh that can be used in openCARP.
//
//  To create the corresponding .msh file, run
//      gmsh -3 bidomain_005cm3.geo
//
//  Converting to CARP directly using meshtool doesn't seem to work, but it can be
//  converted to VTK first from the Gmsh GUI :
//      gmsh bidomain_005cm3.msh
//      Select 'File -> Export', then choose the file extension 'VTK' and save the mesh as 'bidomain_005cm3.vtk'.
//
//  Then, it can be converted to CARP format using
//      meshtool convert -imsh bidomain_005cm3.vtk -ifmt vtk -omsh bidomain_005cm3 -ofmt carp_txt
//  
//  The resulting .elem, .lon, and .pts files is then run with the s1.par file, with
//      openCARP +F s1.par
//
// which gives the error : "L5 : error: Empty stimulus[0] electrode def! Aborting!"
// -----------------------------------------------------------------------------

//Dimensions all in micrometres
lc = 100; //dx = 100um
Lx = 1444;
Ly = 2604;
Lz = 1590;

Point(1) = {0, 0, 0, lc};
Point(2) = {0, 0, Lz, lc};
Point(3) = {Lx, 0, Lz, lc};
Point(4) = {Lx, 0, 0, lc};

Point(5) = {0, Ly, 0, lc};
Point(6) = {0, Ly, Lz, lc};
Point(7) = {Lx, Ly, Lz, lc};
Point(8) = {Lx, Ly, 0, lc};

Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};

Curve Loop(1) = {1, 2, 3, 4};

Plane Surface(1) = {1};

Extrude {0,Ly,0} {Surface{1};}

-------------------------------------------------------------------


and here is the parameters file, s1.par:


# ----------------------------------------
# Output Folder Name
# ----------------------------------------

simID = Test_Output/Bidomain-005cm3-Gmsh

# ----------------------------------------
# Meshname
# ----------------------------------------

meshname = bidomain_005cm3

# ----------------------------------------
# Definition of physics region
# ----------------------------------------
num_phys_regions  = 2
phys_region[0].name = "Intracellular domain"
phys_region[0].ptype = 0
phys_region[0].num_IDs = 1
phys_region[0].ID[0] = 1
phys_region[1].name = "Extracellular domain"
phys_region[1].ptype = 1
phys_region[1].num_IDs = 1
phys_region[1].ID[0] = 1

# ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ---------------
# Definition of ionic model and plugins (IMP) for each unique tagged
# region in the mesh
# ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ------ ----- ----
num_imp_regions = 1
imp_region[0].im = Shannon

num_gregions = 1

# --------------------
# Definition of stimuli
# --------------------

num_stim = 1
#stimulus[0].name = "S1"
#stimulus[0].stimtype = 0
#stimulus[0].strength = 0.0
#stimulus[0].duration = 2.0
#stimulus[0].start    = 0.
#stimulus[0].npls     = 1
#stimulus[0].x0 = 0.0
#stimulus[0].xd = 222.0
#stimulus[0].y0 = 0.0
#stimulus[0].yd = 222.0
#stimulus[0].z0 = 0.0
#stimulus[0].zd = 222.0

stim[0].name = "S1"
stim[0].elec.geom_type = 2
stim[0].elec.p1 = 222.0
stim[0].elec.p0 = 0.0
stim[0].ptcl.duration = 2.0
stim[0].pulse.strength = 300

# -------------------------------
# Definition of general parameters
# -------------------------------
tend = 15.0 # t_end - end of simulation
spacedt = 0.1 # Interval to output data to files
timedt  = 0.1 # Interval to print progress to sceen
vofile = "vm.igb" # Output filename: Prevents invalid read when defined

# --------------------------------------------------
# Definition of solving problem and numerical methods
# --------------------------------------------------
bidomain = 1
parab_solve = 1

---

Hello,

firstly, I would advise you to use the new stimulus definition as follows:

stim[0].elec.p0[0]      = 0.0
stim[0].elec.p1[0]      = 222.0
stim[0].elec.p0[1]      = 0.0
stim[0].elec.p1[1]      = 222.0
stim[0].elec.p0[2]      = 0.0
stim[0].elec.p1[2]      = 222.0

Since I had problems with passing the parameters in curly brackets. I am not completely sure what's happening since in your case the values for the p0 and p1 elements are the same, but this should keep you safe. Secondly, I cannot find the stimulus start in your .par file.

stim[0].ptcl.start      = 0.00

The rest looks good to me. If you continue to have problems, please let me know and I will try to reproduce your results.

---

Thank you for clarifying how to define the stimulus. I tried running it with those values just now, but unfortunately I still got the same error.

If you could try to reproduce the error with my mesh file that would be very helpful. Please let me know if you need anything else beyond what I already sent.

---

I am not that familiar with meshes, but to me your definition looks like a block. You can also define a block using the mesh.Block function provided by carputils and generate it afterwards. Therefore you need to add the following two lines to your run.py script:

geom = mesh.Block(size=(1.444, 2.604, 1.590), resolution=100/1000)
meshname = mesh.generate(geom)

I hope that helps.

---

Yes this would work for the most recent mesh file I posted, but my original issue was that the resolution I hoped to use (2um) was too fine to be generated with mesh.Block or mesher; it resulted in a negative number of tetrahedra and bad array errors when I tried to run it. This is why I switched to using Gmsh. The Gmsh file I posted was just to test out the integration between Gmsh and openCARP, which is what I am currently struggling with (the empty stimulus error).

---

Ok, I got that. I simulated the most recent mesh with the mesh.Block definition:

geom = mesh.Block(size=(1.444, 2.604, 1.590), resolution=args.dx/1000)

For me it is possible with your stimulus definition in a .par file:

stim[0].name            = "S1"
stim[0].elec.geom_type  = 2
stim[0].elec.p0[0]      = 0.0
stim[0].elec.p1[0]      = 222.0
stim[0].elec.p0[1]      = 0.0
stim[0].elec.p1[1]      = 222.0
stim[0].elec.p0[2]      = 0.0
stim[0].elec.p1[2]      = 222.0
stim[0].ptcl.duration   = 2.0
stim[0].pulse.strength  = 300

I do not get an empty stimulus error. I chose several combinations to test, e.g. --dt 10 and --dx 100.

So, this underlines the mesh definition is incorrect. Could you send me the mesh file in opencarp format? Something should be incorrect there. Maybe the coordinates of the electrode and the mesh do not match. Next would be the vtk file and if there are still problems I will install gmsh and check that.

---

Hi Moritz, thank you for your reply. In your testing did you try --dx 2 yet? If not, could you please try this value? The mesh and stimulus do work for all dx values greater than this, but I am interested in a resolution of 0.002 mm.

---

However, since I can't get any resolution working from an external mesh, I guess either way it makes sense to share the opencarp mesh files with you. Where can I send it?

---

Sorry for the delayed response. The problem is within the conversion of gmsh to openCARP readable mesh.

There are two ways which fixed the issue for me you could try:
1) you open the file with gmsh and add a volume "Physical groups" -> "Geometry" -> "Physical groups" -> "Add" -> "Volume" by clicking on the volume before you safe the mesh as vtk-file.

2) use a short bash script:
GMSH="/opt/homebrew/bin/gmsh"
$GMSH -parse_and_exit bidomain_005cm3.geo
sed -i 's/CellEntityIds/elemTags/g' bidomain_005cm3.vtk
meshtool convert -imsh=bidomain_005cm3.vtk -omsh=bidomain_005cm3 -ofmt=carp_txt

The mesh name is given by the path to your mesh files:
meshname = 'meshes/test/bidomain_005cm3'

After the conversion of the gmsh file you need to adjust your run.py and set the:
IntraTags = [27] and ExtraTags = [27]. This is the default, which is set if you add a volume the gmsh geometry.

I hope this solves the problem.

---

Hi Moritz, thank you for the response. After following those two steps, the test mesh (dx 100um) is working with openCARP now!

I am unfortunately having issues now with generating my original mesh in Gmsh, which has a dx of 2um.  At the start of this thread, Axel had suggested I use Gmsh because Mesher could not handle large meshes, however, despite using Gmsh, the process keeps getting killed. I am not sure if you can help with this as Gmsh is outside of your domain. If you happen to have any suggestions though, please let me know!

Either way, thanks again for your help.