Add region mesh style defined by a closed STL or surf file surface, 2d and 3d

cmake/common/set/sparta_cmake_defaults.cmake

@@ -56,6 +56,7 @@ if(SPARTA_ENABLE_TESTING)
       "explicit2implicit"
       "mfp_mct"
       "optmove"
+      "region_mesh"
       "torque")
 
   set(SPARTA_DISABLED_TESTS

data/README

@@ -19,11 +19,15 @@ Files that end in ".tce" list reactions gas-phase chemistry for the
 Total Collision Energy (TCE) model and are read by the react command.
 
 Files that start with "sdata." contain lists of surface elements and
-are read by the read_surf command.
+are read by the read_surf command.  They can also be read by the region
+command with its mesh style, which defines the region as the volume
+enclosed by the surface.
 
 Files that end in ".stl" are STL-formatted surface files which can be
 converted to SPARTA formatted surface files via the stl2surf.py tool
-in the tools directory.
+in the tools directory.  They can also be read directly by the region
+command with its mesh style, which defines the region as the volume
+enclosed by the surface.
 
 Files that end in ".surf" list reactions for surface chemistry
 and are read by the surf_react command.
@@ -35,6 +39,6 @@ species = doc/species.html
 species.vib = doc/species.html
 vss/vhs/hs = doc/collide.html
 tce = doc/react.html
-sdata = doc/read_surf.html
-stl = doc/Section_tools.html#stl2surf
+sdata = doc/read_surf.html and doc/region.html
+stl = doc/Section_tools.html#stl2surf and doc/region.html
 surf = doc/surf_react.html

doc/region.html

@@ -17,7 +17,7 @@ <H3>region command
 </PRE>
 <UL><LI>ID = user-assigned name for the region 
 
-<LI>style = <I>block</I> or <I>cylinder</I> or <I>plane</I> or <I>sphere</I> or <I>union</I> or <I>intersect</I> or <I>block/kk</I> or <I>cylinder/kk</I> or <I>plane/kk</I> or <I>sphere/kk</I> or <I>union/kk</I> or <I>intersect/kk</I> 
+<LI>style = <I>block</I> or <I>cylinder</I> or <I>mesh</I> or <I>plane</I> or <I>sphere</I> or <I>union</I> or <I>intersect</I> or <I>block/kk</I> or <I>cylinder/kk</I> or <I>plane/kk</I> or <I>sphere/kk</I> or <I>union/kk</I> or <I>intersect/kk</I> 
 
 <PRE>  <I>block</I> or <I>block/kk</I> args = xlo xhi ylo yhi zlo zhi
     xlo,xhi,ylo,yhi,zlo,zhi = bounds of block in all dimensions (distance units)
@@ -26,6 +26,9 @@ <H3>region command
     c1,c2 = coords of cylinder axis in other 2 dimensions (distance units)
     radius = cylinder radius (distance units)
     lo,hi = bounds of cylinder in dim (distance units)
+  <I>mesh</I> args = file
+    file = name of file with a closed surface, in STL format or in the
+           format the <A HREF = "read_surf.html">read_surf</A> command reads
   <I>plane</I> or <I>plane/kk</I> args = px py pz nx ny nz
     px,py,pz = point on the plane (distance units)
     nx,ny,nz = direction normal to plane (distance units)
@@ -41,11 +44,20 @@ <H3>region command
 </PRE>
 <LI>zero or more keyword/value pairs may be appended 
 
-<LI>keyword = <I>side</I> 
+<LI>keyword = <I>side</I> or <I>origin</I> or <I>trans</I> or <I>scale</I> or <I>rotate</I> 
 
 <PRE>  <I>side</I> value = <I>in</I> or <I>out</I>
     <I>in</I> = the region is inside the specified geometry
-    <I>out</I> = the region is outside the specified geometry 
+    <I>out</I> = the region is outside the specified geometry
+  <I>origin</I> values = Ox Oy Oz (only for style <I>mesh</I>)
+    Ox,Oy,Oz = point that <I>scale</I> and <I>rotate</I> operate around (distance units)
+  <I>trans</I> values = Dx Dy Dz (only for style <I>mesh</I>)
+    Dx,Dy,Dz = distance to translate the mesh by (distance units)
+  <I>scale</I> values = Sx Sy Sz (only for style <I>mesh</I>)
+    Sx,Sy,Sz = scale factor to apply to the mesh in each dimension
+  <I>rotate</I> values = theta Rx Ry Rz (only for style <I>mesh</I>)
+    theta = angle to rotate the mesh by (degrees)
+    Rx,Ry,Rz = axis of rotation 
 </PRE>
 
 </UL>
@@ -54,7 +66,10 @@ <H3>region command
 <PRE>region 1 block -3.0 5.0 INF 10.0 INF INF
 region 2 sphere 0.0 0.0 0.0 5 side out
 region void cylinder y 2 3 5 -5.0 INF
-region outside union 4 side1 side2 side3 side4 
+region outside union 4 side1 side2 side3 side4
+region flask mesh flask.stl
+region duct mesh duct.stl scale 0.001 0.001 0.001 side out
+region shell mesh data.circle side out 
 </PRE>
 <P><B>Description:</B>
 </P>
@@ -82,6 +97,53 @@ <H3>region command
 y-direction located at x = 2.0 and z = 3.0, with a radius of 5.0, and
 extending in the y-direction from -5.0 to infinity.
 </P>
+<P>The <I>mesh</I> style defines the region as the volume enclosed by a closed
+surface which is read from a file.  Two file formats are supported and
+SPARTA detects which one a file is from its contents, not from its name:
+</P>
+<P>STL (stereolithography) files, in either of the two variants of that
+format, plain text ("ASCII") or binary.  An STL file always stores a 3d
+triangulated surface.
+</P>
+<P>Surface files in the format the <A HREF = "read_surf.html">read_surf</A> command
+reads, which store a 3d surface as triangles or a 2d surface as line
+segments.  As for read_surf, a <I>Points</I> section is optional, and a
+per-surface type column is skipped if it is present.  Per-surface custom
+attributes are not supported.
+</P>
+<P>Only the geometry in the file is used.  Surface types, surface IDs, and
+the direction of the normals are all ignored.  Every processor reads and
+stores the entire mesh, so this style is best suited to small and medium
+sized meshes.
+</P>
+<P>A point is tested against a 3d mesh by casting a ray from the point in
+the +x direction and counting how many triangles it crosses.  An odd
+count means the point is inside the surface.  A 2d mesh is a closed
+curve in the x-y plane, so the ray is cast in that plane and crossings
+of the line segments are counted instead.  The z coordinate of the point
+does not matter for a 2d mesh, which means the region it defines is the
+curve extruded to infinity in both z directions.  A 2d mesh can thus be
+used in a 3d simulation as well, where it defines a prism.
+</P>
+<P>Either way this requires the surface to be closed, and SPARTA prints a
+warning if the elements in the file do not enclose a volume (or an area
+in 2d), in which case the inside/outside test gives arbitrary answers.
+As for the other styles, a point that lies exactly on the surface is
+considered to be inside it, and the <I>side</I> keyword selects whether the
+region is the volume inside or outside the mesh.
+</P>
+<P>The <I>origin</I>, <I>trans</I>, <I>scale</I>, and <I>rotate</I> keywords can only be used
+with the <I>mesh</I> style.  They transform the mesh after it is read, which
+is useful when the file uses different units or is positioned
+differently than the simulation box.  They are applied in the order they
+appear in the command, and follow the same conventions as the
+<A HREF = "read_surf.html">read_surf</A> keywords of the same name.  <I>Trans</I>
+translates the mesh by the specified distance.  <I>Scale</I> multiplies the
+distance of every vertex from the <I>origin</I> point by the specified
+factors.  <I>Rotate</I> rotates the mesh by theta degrees around the axis
+Rx,Ry,Rz which passes through the <I>origin</I> point.  The <I>origin</I> point is
+0,0,0 by default and is moved along by any preceding <I>trans</I>.
+</P>
 <P>The <I>union</I> style creates a region consisting of the volume of all the
 listed regions combined.  The <I>intersect</I> style creates a region
 consisting of the volume that is common to all the listed regions.
@@ -127,11 +189,28 @@ <H3>region command
 </P>
 <HR>
 
-<P><B>Restrictions:</B> none
+<P><B>Restrictions:</B>
+</P>
+<P>An STL file always describes a 3d surface, so a 2d region can only be
+defined with a <A HREF = "read_surf.html">read_surf</A> format file containing lines.
+Unlike read_surf, the <I>mesh</I> style does not require the dimensionality
+of the file to match the <A HREF = "dimension.html">dimension</A> of the simulation,
+since as noted above regions in SPARTA are 3d objects in either case.
+</P>
+<P>There is no <I>mesh/kk</I> style.  The
+<A HREF = "create_particles.html">create_particles</A> command tests points against a
+region on the host, so a <I>mesh</I> region can be used with
+create_particles/kk, either on its own or as one of the regions listed
+by a <I>union</I> or <I>intersect</I> region.  Commands that test regions on the
+device, such as <A HREF = "fix_emit_face.html">fix emit/face</A> and <A HREF = "fix_ave_histo.html">fix ave/histo</A>,
+will generate an error if they are given a <I>mesh</I> region when running
+with the KOKKOS package, whether it is given to them directly or nested
+inside a <I>union</I> or <I>intersect</I> region.
 </P>
 <P><B>Related commands:</B>
 </P>
-<P><A HREF = "dump_modify.html">dump_modify</A>
+<P><A HREF = "create_particles.html">create_particles</A>, <A HREF = "read_surf.html">read_surf</A>,
+<A HREF = "dump_modify.html">dump_modify</A>
 </P>
 <P><B>Default:</B>
 </P>

doc/region.txt

@@ -13,14 +13,17 @@ region command :h3
 region ID style args keyword value ... :pre
 
 ID = user-assigned name for the region :ulb,l
-style = {block} or {cylinder} or {plane} or {sphere} or {union} or {intersect} or {block/kk} or {cylinder/kk} or {plane/kk} or {sphere/kk} or {union/kk} or {intersect/kk} :l
+style = {block} or {cylinder} or {mesh} or {plane} or {sphere} or {union} or {intersect} or {block/kk} or {cylinder/kk} or {plane/kk} or {sphere/kk} or {union/kk} or {intersect/kk} :l
   {block} or {block/kk} args = xlo xhi ylo yhi zlo zhi
     xlo,xhi,ylo,yhi,zlo,zhi = bounds of block in all dimensions (distance units)
   {cylinder} or {cylinder/kk} args = dim c1 c2 radius lo hi
     dim = {x} or {y} or {z} = axis of cylinder
     c1,c2 = coords of cylinder axis in other 2 dimensions (distance units)
     radius = cylinder radius (distance units)
     lo,hi = bounds of cylinder in dim (distance units)
+  {mesh} args = file
+    file = name of file with a closed surface, in STL format or in the
+           format the "read_surf"_read_surf.html command reads
   {plane} or {plane/kk} args = px py pz nx ny nz
     px,py,pz = point on the plane (distance units)
     nx,ny,nz = direction normal to plane (distance units)
@@ -34,18 +37,30 @@ style = {block} or {cylinder} or {plane} or {sphere} or {union} or {intersect} o
     N = # of regions to follow, must be 2 or greater
     reg-ID1,reg-ID2, ... = IDs of regions to intersect :pre
 zero or more keyword/value pairs may be appended :l
-keyword = {side} :l
+keyword = {side} or {origin} or {trans} or {scale} or {rotate} :l
   {side} value = {in} or {out}
     {in} = the region is inside the specified geometry
-    {out} = the region is outside the specified geometry :pre
+    {out} = the region is outside the specified geometry
+  {origin} values = Ox Oy Oz (only for style {mesh})
+    Ox,Oy,Oz = point that {scale} and {rotate} operate around (distance units)
+  {trans} values = Dx Dy Dz (only for style {mesh})
+    Dx,Dy,Dz = distance to translate the mesh by (distance units)
+  {scale} values = Sx Sy Sz (only for style {mesh})
+    Sx,Sy,Sz = scale factor to apply to the mesh in each dimension
+  {rotate} values = theta Rx Ry Rz (only for style {mesh})
+    theta = angle to rotate the mesh by (degrees)
+    Rx,Ry,Rz = axis of rotation :pre
 :ule
 
 [Examples:]
 
 region 1 block -3.0 5.0 INF 10.0 INF INF
 region 2 sphere 0.0 0.0 0.0 5 side out
 region void cylinder y 2 3 5 -5.0 INF
-region outside union 4 side1 side2 side3 side4 :pre
+region outside union 4 side1 side2 side3 side4
+region flask mesh flask.stl
+region duct mesh duct.stl scale 0.001 0.001 0.001 side out
+region shell mesh data.circle side out :pre
 
 [Description:]
 
@@ -73,6 +88,53 @@ third example above specifies a cylinder with its axis in the
 y-direction located at x = 2.0 and z = 3.0, with a radius of 5.0, and
 extending in the y-direction from -5.0 to infinity.
 
+The {mesh} style defines the region as the volume enclosed by a closed
+surface which is read from a file.  Two file formats are supported and
+SPARTA detects which one a file is from its contents, not from its name:
+
+STL (stereolithography) files, in either of the two variants of that
+format, plain text ("ASCII") or binary.  An STL file always stores a 3d
+triangulated surface.
+
+Surface files in the format the "read_surf"_read_surf.html command
+reads, which store a 3d surface as triangles or a 2d surface as line
+segments.  As for read_surf, a {Points} section is optional, and a
+per-surface type column is skipped if it is present.  Per-surface custom
+attributes are not supported.
+
+Only the geometry in the file is used.  Surface types, surface IDs, and
+the direction of the normals are all ignored.  Every processor reads and
+stores the entire mesh, so this style is best suited to small and medium
+sized meshes.
+
+A point is tested against a 3d mesh by casting a ray from the point in
+the +x direction and counting how many triangles it crosses.  An odd
+count means the point is inside the surface.  A 2d mesh is a closed
+curve in the x-y plane, so the ray is cast in that plane and crossings
+of the line segments are counted instead.  The z coordinate of the point
+does not matter for a 2d mesh, which means the region it defines is the
+curve extruded to infinity in both z directions.  A 2d mesh can thus be
+used in a 3d simulation as well, where it defines a prism.
+
+Either way this requires the surface to be closed, and SPARTA prints a
+warning if the elements in the file do not enclose a volume (or an area
+in 2d), in which case the inside/outside test gives arbitrary answers.
+As for the other styles, a point that lies exactly on the surface is
+considered to be inside it, and the {side} keyword selects whether the
+region is the volume inside or outside the mesh.
+
+The {origin}, {trans}, {scale}, and {rotate} keywords can only be used
+with the {mesh} style.  They transform the mesh after it is read, which
+is useful when the file uses different units or is positioned
+differently than the simulation box.  They are applied in the order they
+appear in the command, and follow the same conventions as the
+"read_surf"_read_surf.html keywords of the same name.  {Trans}
+translates the mesh by the specified distance.  {Scale} multiplies the
+distance of every vertex from the {origin} point by the specified
+factors.  {Rotate} rotates the mesh by theta degrees around the axis
+Rx,Ry,Rz which passes through the {origin} point.  The {origin} point is
+0,0,0 by default and is moved along by any preceding {trans}.
+
 The {union} style creates a region consisting of the volume of all the
 listed regions combined.  The {intersect} style creates a region
 consisting of the volume that is common to all the listed regions.
@@ -118,10 +180,28 @@ effectively.
 
 :line
 
-[Restrictions:] none
+[Restrictions:]
+
+An STL file always describes a 3d surface, so a 2d region can only be
+defined with a "read_surf"_read_surf.html format file containing lines.
+Unlike read_surf, the {mesh} style does not require the dimensionality
+of the file to match the "dimension"_dimension.html of the simulation,
+since as noted above regions in SPARTA are 3d objects in either case.
+
+There is no {mesh/kk} style.  The
+"create_particles"_create_particles.html command tests points against a
+region on the host, so a {mesh} region can be used with
+create_particles/kk, either on its own or as one of the regions listed
+by a {union} or {intersect} region.  Commands that test regions on the
+device, such as "fix emit/face"_fix_emit_face.html and "fix
+ave/histo"_fix_ave_histo.html, will generate an error if they are given
+a {mesh} region when running with the KOKKOS package, whether it is
+given to them directly or nested inside a {union} or {intersect}
+region.
 
 [Related commands:]
 
+"create_particles"_create_particles.html, "read_surf"_read_surf.html,
 "dump_modify"_dump_modify.html
 
 [Default:]

examples/region_mesh/air.species

@@ -0,0 +1,24 @@
+# Species data
+
+# ID
+# Molwt (amu)
+# Molmass (kg)
+# Rotational dof
+# RotRel
+# Vibrational dof
+# VibRel
+# VibTemp (K)
+# species wt
+# charge
+
+O2  32.00    5.31E-26  2    0.2   2    5.58659E-5    2256.0    1.0      0.0
+N2  28.016   4.65E-26  2    0.2   2    1.90114E-5    3371.0    1.0      0.0
+O   16.00    2.65E-26  0    0.0   0    0.0           0.0       1.0      0.0
+N   14.008  2.325E-26  0    0.0   0    0.0           0.0       1.0      0.0
+NO  30.008   4.98E-26  2    0.2   2    7.14285E-4    2719.0    1.0      0.0
+O2+ 32.00    5.31E-26  2    0.2   2    5.58659E-5    2256.0    1.0      1.0
+N2+ 28.016   4.65E-26  2    0.2   2    1.90114E-5    3371.0    1.0      1.0
+O+  16.00    2.65E-26  0    0.0   0    0.0           0.0       1.0      1.0
+N+  14.008  2.325E-26  0    0.0   0    0.0           0.0       1.0      1.0
+NO+ 30.008   4.98E-26  2    0.2   2    7.14285E-4    2719.0    1.0      1.0
+e   0.001  9.10938188E-31 0 0.0   0    0.0           0.0       1.0     -1.0

examples/region_mesh/air.vss

@@ -0,0 +1,22 @@
+# VSS collision model parameters for each species
+
+# diameter (mass)
+# omega (unitless)
+# tref (temperature)
+# alpha (unitless)
+# Zrotinf (unitless)
+# T* (temperature)
+# C1 (temperature)
+# C2 (temperature^(1/3))
+
+O2   3.96E-10    0.77  273.15  1.4   16.5 113.5 56.5 153.5
+N2   4.07E-10    0.74  273.15  1.6   18.1 91.5 9.1 220.0
+O    3.0E-10     0.80  273.15  1.0   0.0 0.0 0.0 0.0
+N    3.0E-10     0.80  273.15  1.0   0.0 0.0 0.0 0.0
+NO   4.0E-10     0.80  273.15  1.0   7.5 119.0 9.10 220.00
+O2+  3.96E-10    0.77  273.15  1.4   16.5 113.5 56.5 153.5
+N2+  4.07E-10    0.74  273.15  1.6   18.1 91.5 9.1 220.0
+O+   3.0E-10     0.80  273.15  1.0   0.0 0.0 0.0 0.0
+N+   3.0E-10     0.80  273.15  1.0   0.0 0.0 0.0 0.0
+NO+  4.0E-10     0.80  273.15  1.0   7.5 119.0 9.10 220.00
+e    7.0E-13     0.50  273.15  1.0   0.0 0.0 0.0 0.0