diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index c1e730a37e0..c810d5b0334 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -83,6 +83,7 @@ OPT. * :doc:`gjf ` * :doc:`gld ` * :doc:`gle ` + * :doc:`graphics ` * :doc:`gravity (ko) ` * :doc:`grem ` * :doc:`halt ` diff --git a/doc/src/JPG/fix-graphics-example.png b/doc/src/JPG/fix-graphics-example.png new file mode 100644 index 00000000000..d9b1d4a5cbf Binary files /dev/null and b/doc/src/JPG/fix-graphics-example.png differ diff --git a/doc/src/Modify_fix.rst b/doc/src/Modify_fix.rst index 18b51d6e137..fd089efa5db 100644 --- a/doc/src/Modify_fix.rst +++ b/doc/src/Modify_fix.rst @@ -129,6 +129,8 @@ derived class. See ``src/fix.h`` for additional details. +---------------------------+--------------------------------------------------------------------------------------------+ | thermo | compute quantities for thermodynamic output (optional) | +---------------------------+--------------------------------------------------------------------------------------------+ +| image | pass lists of graphics objects and their parameters to :doc:`dump image fix ` | ++---------------------------+--------------------------------------------------------------------------------------------+ Typically, only a small fraction of these methods are defined for a particular fix. Setmask is mandatory, as it determines when the fix diff --git a/doc/src/dump_image.rst b/doc/src/dump_image.rst index ec769b5e8f5..ed77daef314 100644 --- a/doc/src/dump_image.rst +++ b/doc/src/dump_image.rst @@ -81,12 +81,15 @@ Syntax axes = *yes* or *no* = do or do not draw xyz axes lines next to simulation box length = length of axes lines as fraction of respective box lengths diam = diameter of axes lines as fraction of shortest box length - *region* values = region-ID color drawstyle [npoints (optional) diameter (optional)] + *region* values = region-ID color drawstyle [opacity (optional) npoints (optional) diameter (optional)] region-ID = ID of the region to render color = color name for region graphics - drawstyle = *filled* or *frame* or *points* + drawstyle = *filled* or *transparent* or *frame* or *points* *filled* = render region as a filled object, with optional open faces + *transparent* = same as *filled* but has selectable opacity *frame* = render region as a wireframe (like box or subbox) + *points* = fill region with spheres at random locations + opacity = level of opacity (from 0.0 to 1.0, only for drawstyle *transparent*) npoints = number of attempted points (only for drawstyle *points*) diameter = diameter of wireframe or points (only for drawstyles *frame* and *points*) *subbox* values = lines diam = draw outline of processor subdomains @@ -114,7 +117,7 @@ Syntax dump_modify dump-ID keyword values ... * these keywords apply only to the *image* and *movie* styles and are documented on this page -* keyword = *acolor* or *adiam* or *amap* or *gmap* or *backcolor* or *bcolor* or *bdiam* or *bitrate* or *boxcolor* or *color* or *framerate* or *gmap* +* keyword = *acolor* or *adiam* or *amap* or *gmap* or *atrans* or *backcolor* or *bcolor* or *bdiam* or *btrans* or *bitrate* or *boxcolor* or *color* or *framerate* or *axestrans* or *boxtrans* or *subboxtrans* * see the :doc:`dump modify ` doc page for more general keywords .. parsed-literal:: @@ -145,6 +148,9 @@ Syntax color = name of color used for that subset of values entry = color (for sequential style) color = name of color used for a bin of values + *atrans* args = type transparency + type = atom type (numeric or type label) or range of numeric types (see below) + transparency = transparency of atoms of that type (value between 0 (invisible) and 1 (fully opaque)) *backcolor* arg = color color = name of color for background *bcolor* args = type color @@ -153,13 +159,25 @@ Syntax *bdiam* args = type diam type = bond type (numeric or type label) or range of numeric types (see below) diam = diameter of bonds of that type (distance units) - *bitrate* arg = rate - rate = target bitrate for movie in kbps + *btrans* args = type transparency + type = bond type (numeric or type label) or range of numeric types (see below) + transparency = transparency of bonds of that type (value between 0 (invisible) and 1 (fully opaque)) + *axestrans* arg = transparency + transparency = transparency for axes lines (value between 0 (invisible) and 1 (fully opaque)) + *boxtrans* arg = transparency + transparency = transparency for simulation box lines (value between 0 (invisible) and 1 (fully opaque)) *boxcolor* arg = color color = name of color for simulation box lines and processor subdomain lines + *subboxtrans* arg = transparency + transparency = transparency for simulation subbox lines (value between 0 (invisible) and 1 (fully opaque)) *color* args = name R G B name = name of color R,G,B = red/green/blue numeric values from 0.0 to 1.0 + *fcolor* args = fix-ID color + fix-ID = ID of the fix + color = name of color for image objects provided by this fix + *bitrate* arg = rate + rate = target bitrate for movie in kbps *framerate* arg = fps fps = frames per second for movie *gmap* args = identical to *amap* args @@ -376,6 +394,9 @@ the mass of both atoms of a pair within the bond cutoff is lower than 3 atomic mass units, a bond is **not** drawn; this prohibits displaying unwanted hydrogen-hydrogen bonds for alkyl or alcohol groups or for water with typical cutoffs suitable for displaying covalent bonds. +For ReaxFF it is also possible to visualize bonds as they are computed +through using :doc:`fix reaxff/bonds ` with the +*fix* keyword (see below). ---------- @@ -504,46 +525,77 @@ change this via the dump_modify command. ---------- -The *fix* keyword can be used with a :doc:`fix ` that produces -objects to be drawn. - -The *fflag1* and *fflag2* settings are numerical values which are -passed to the fix to affect how the drawing of its objects is done. -See the individual fix page for a description of what these -parameters mean for a particular fix. - -The only setting currently allowed for the *color* value is *type*, -which will color the fix objects according to their type. By default -the mapping of types to colors is as follows: +.. versionchanged:: TBD -* type 1 = red -* type 2 = green -* type 3 = blue -* type 4 = yellow -* type 5 = aqua -* type 6 = cyan + Support for several fix styles added and more flexible color selection -and repeats itself for types > 6. There is not yet an option to -change this via the dump_modify command. +The *fix* keyword can be used with a :doc:`fix ` that produces +objects to be drawn. Below is a list of supported fixes: + +* :doc:`fix graphics ` +* :doc:`fix indent ` +* :doc:`fix smd/wall_surface ` +* :doc:`fix wall/lj93 ` +* :doc:`fix wall/lj126 ` +* :doc:`fix wall/lj1043 ` +* :doc:`fix wall/colloid ` +* :doc:`fix wall/gran ` +* :doc:`fix wall/harmonic ` +* :doc:`fix wall/harmonic/outside ` +* :doc:`fix wall/lepton ` +* :doc:`fix wall/morse ` +* :doc:`fix wall/reflect ` +* :doc:`fix wall/reflect/stochastic ` +* :doc:`fix wall/table ` + +The fix keyword may be used multiple times to include visualizations of +graphics objects from multiple fixes. The fix keyword is followed by +the :doc:`fix ID ` of the fix, the color style setting and two +numerical values *fflag1* and *fflag2*. + +The color style may be either *type*, *element*, or *const*. The first +two will use the same color as assigned to the corresponding atom type +and thus it depends on the fix which atom type it associates with any +object. Often this will be atom type 1. For the *const* type a +constant color will be used that can be changed with a *dump_modify +fcolor* command (see below). By default the constant color will be +"red" (same as the default color for atom type 1). + +The *fflag1* and *fflag2* settings are numerical values which are used +by *dump image* to adjust how the drawing of the objects communicated by +the fix is done. See the documentation of the individual fixes for a +description of what these parameters mean for the graphics objects +provided by those fixes. ---------- .. versionadded:: 10Sep2025 +.. versionchanged:: TBD + + style *transparency* was added + The *region* keyword can be used to create a graphical representation of a :doc:`region `. This can be helpful in debugging the location and extent of regions, especially when those have parameters controlled by variables. Three styles of representing a region are available: -*filled*, *frame*, and *points*. With style *filled* the surface of the -region is drawn. For region styles that support open faces, surfaces -are not drawn for such open faces. Draw style *frame* represents the -region with a mesh of "wires". The diameter of these "wires" can be -set. Unlike with the *filled* style, you can see what is *inside* the -region with this draw style. The third draw style *points* generates a -random point cloud inside the simulation box and draws only those points -that are within the region. Draw styles *filled* and *frame* support -only "primitive" region style (no unions or intersections), but the -*points* draw style supports all region styles. +*filled*\, *transparency*\, *frame*\, and *points*. With style *filled* +the surface of the region is triangulated and drawn. For region styles +that support open faces, surfaces for such open faces are skipped. The +style *transparent* is like *filled* but takes an additional parameter +in the range of 0.0 to 1.0 that defines the opacity and thus allows to +see what is inside the region. Draw style *frame* represents the region +with a mesh of "wires". The diameter of these "wires" can be set. +Unlike with the *filled* style and similar to the *transparent* style, +you can see what is *inside* the region with this draw style. The third +draw style, *points*\, generates a random point cloud inside the +simulation box and draws only those points that are within the region. +Draw styles *filled*\, *transparent*\, and *frame* support only +"primitive" region styles (no unions or intersections), but the *points* +draw style supports *all* region styles. + +Recommended transparency settings are the values of 0.25, 0.5, or 0.75 +when used in combination with *fsaa on*. ---------- @@ -1029,6 +1081,35 @@ pre-defined color names with new RBG values. ---------- +**Transparency settings for atoms bonds and standard visualization objects** + +.. versionadded:: TBD + +Various graphical objects in *dump image* output can be rendered in a +transparent fashion using the so-called screen-door transparency method. +This means that only a subset of pixels for a graphical object are +written to the image. This can be controlled with various +*dump\_modify* settings: *atrans* for atoms, *btrans* for bonds, +*axestrans* for axes lines, *boxtrans* for the simulation box, and +*subboxtrans* for the subdomain box lines. The transparency value +must be between 0.0 (invisible) and 1.0 (fully opaque). The default +setting for all is 1.0. + +Recommended transparency settings are the values of 0.25, 0.5, or 0.75 +when used in combination with *fsaa on*. + +---------- + +.. versionadded:: TBD + +The *fcolor* keyword sets the color of any image objects created by a +fix. The first argument is the fix ID used with the *dump image fix* +command and the second argument is the color name. The color name can +be any of the 140 pre-defined colors (see below) or a color name defined +by the *dump_modify color* option. + +---------- + The *framerate* keyword can be used with the :doc:`dump movie ` command to define the duration of the resulting movie file. Movie files written by the dump *movie* command have a default @@ -1122,6 +1203,7 @@ The defaults for the dump image and dump movie keywords are as follows: * subbox no 0.0 * shiny = 1.0 * ssao = no +* fsaa = no ---------- @@ -1130,14 +1212,18 @@ The defaults for the dump_modify keywords specific to dump image and dump movie * acolor = \* red/green/blue/yellow/aqua/cyan * adiam = \* 1.0 * amap = min max cf 0.0 2 min blue max red +* atrans = 1.0 * backcolor = black * bcolor = \* red/green/blue/yellow/aqua/cyan * bdiam = \* 0.5 -* bitrate = 2000 +* btrans = 1.0 * boxcolor = yellow +* axestrans = 1.0 +* boxtrans = 1.0 +* subboxtrans = 1.0 * color = 140 color names are pre-defined as listed below +* bitrate = 2000 * framerate = 24 -* fsaa = no * gmap = min max cf 0.0 2 min blue max red ---------- diff --git a/doc/src/fix.rst b/doc/src/fix.rst index cda13dcf7c3..2a71968ba1e 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -262,6 +262,7 @@ accelerated styles exist. * :doc:`gjf ` - statistically correct Langevin temperature control using the GJ methods * :doc:`gld ` - generalized Langevin dynamics integrator * :doc:`gle ` - generalized Langevin equation thermostat +* :doc:`graphics ` - add graphics objects to :doc:`dump image ` output * :doc:`gravity ` - add gravity to atoms in a granular simulation * :doc:`grem ` - implements the generalized replica exchange method * :doc:`halt ` - terminate a dynamics run or minimization diff --git a/doc/src/fix_graphics.rst b/doc/src/fix_graphics.rst new file mode 100644 index 00000000000..b3247e2f976 --- /dev/null +++ b/doc/src/fix_graphics.rst @@ -0,0 +1,230 @@ +.. index:: fix graphics + +fix graphics command +==================== + +Syntax +"""""" + +.. code-block:: LAMMPS + + fix ID group-ID graphics Nevery keyword args ... + +* ID, group-ID are documented in :doc:`fix ` command +* graphics = style name of this fix command +* Nevery = update graphics information every this many time steps +* one or more keyword/args pairs may be appended +* keyword = *sphere* or *cylinder* or *arrow* or *progbar* + + .. parsed-literal:: + + *sphere* args = type x y z R + type = an atom type value to select the color of the sphere + x, y, z = position of the center of the sphere (distance units) + R = sphere radius (distance units) + any of x, y, z, and R can be a variable (see below) + *cylinder* args = type x1 y1 z1 x2 y2 z2 R + type = an atom type value to select the color of the sphere + x1, y1, z1, x2, y2, z2 = positions of the centers at the two ends of the cylinder (distance units) + R = cylinder radius (distance units) + any of x1, y1, z1, x2, y2, z2, and R can be a variable (see below) + *arrow* args = type x1 y1 z1 x2 y2 z2 R ratio + type = an atom type value to select the color of the sphere + x1, y1, z1, x2, y2, z2 = positions of the centers at the tip and the bottom of the arrow (distance units) + R = cylinder radius (distance units) + ratio = tip to body ratio (unitless) + any of x1, y1, z1, x2, y2, z2, and R can be a variable (see below) + *progbar* args = type1 type2 dim x y z length R ratio tics + type1 = an atom type value to select the color of the progress bar body and the tics + type2 = an atom type value to select the color of the progress indicator + dim = *x* or *y* or *z*, direction of the progress bar + x, y, z = position of the progress bar center (distance units) + length = length of progress bar (distance units) + R = cylinder radius (distance units) + ratio = progress status (unitless) + tics = number of tics (unitless) + only the progress ratio value can be a variable (see below) + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix 1 all graphics 100 sphere 1 0.0 0.0 15.0 3.0 sphere 2 0.0 0.0 5.0 1.0 + fix 1 all graphics 1000 sphere 1 v_x v_y 0.0 v_radius cylinder 1 v_x v_y 0.0 v_x v_y 10.0 3.0 + fix 2 all graphics 100 progbar 3 1 z 0.012 -0.012 0.0025 0.03 0.0003 v_prog 10 + +Description +""""""""""" + +.. versionadded:: TBD + +This fix allows to add arbitrary objects to images rendered with +:doc:`dump image ` using the *fix* keyword. + +The *Nevery* keyword determines how often the graphics object data is +updated. This should be the same value as the corresponding *N* +parameter of the :doc:`dump ` image command. LAMMPS will stop +with an error message if the settings for this fix and the dump command +are not compatible. + +Available graphics objects are (see above for exact command line syntax): + +- *sphere* - a sphere defined by its center location and its radius +- *cylinder* - a cylinder defined by its two center endpoints and its radius +- *arrow* - a cylinder with a cone at one side (see note below) +- *progbar* - progress bar a long a selected axis and with optional tick marks + +The *type* quantity determines the color of the object. Its represents +an *atom* type and the object will be colored the same as the +corresponding atom type when the *type* coloring scheme is used in the +:doc:`dump image fix ` command is used. The color may also +be that of the atom type's element or just a globally set constant color +for *all* objects of this fix instance, which can be changed using a +:doc:`dump modify fcolor ` command. For the *progbar* +object *two* atom type values must be specified. + +The *x*\, *y*\, and *z* parameters correspond to the position of the +center of the object (*sphere* and *progbar*). *x1*\, *y1*\, and *z1* as +well as *x2*\, *y2*\, and *z2* are instead representing the top and +bottom position of a graphics object (*cylinder* and *arrow*). The *R* +parameter determines the radius. + +The *progbar* object has four additional parameters: *dim* sets the +direction of the progress bar, "x", "y", or "z"; *length* sets the +length of the entire object; *ratio* sets the ratio of progress and +is expected to be between 0.0 and 1.0 (larger or smaller values will +be reset to 1.0 or 0.0, respectively); and *tics* determines the number +of tics shown on the progress bar, this must be a number between 0 and 20. +Unlike for the other graphics objects, all settings except for *ratio* +are fixed and cannot be a variable reference. + +.. admonition:: Work in progress notice + :class: note + + The *arrow* object is currently composed of two cylinders since the + :doc:`dump image ` render implementation is missing a + primitive to render a cone. This fix will be updated when that + functionality becomes available. + +---------------------- + +Many of the quantities defining a graphics object can be specified as an +equal-style :doc:`variable `, namely *x*, *y*, *z*, or *R* for +a *sphere* or namely *x1*, *y1*, *z1*, *x2*, *y2*, *z2*, or *R* for a +*cylinder*. If any of these values is a variable, it should be +specified as `v_name`, where `name` is the variable name. In this case, +the variable will be evaluated each *Nevery* timestep, and its value +used to define the indenter geometry. + +Note that equal-style variables can specify formulas with various +mathematical functions, and include :doc:`thermo_style ` +command keywords for the simulation box parameters and timestep and +elapsed time. Thus it is easy to specify graphics object properties +like position, orientation, radius or more that change as a function of +time or span consecutive runs in a continuous fashion. For the latter, +see the *start* and *stop* keywords of the :doc:`run ` command and +the *elaplong* keyword of :doc:`thermo_style custom ` for +details. + +For example, if a sphere's x-position is specified as v_x, then this +variable definition will keep its center at a relative position in the +simulation box, 1/4 of the way from the left edge to the right edge, +even if the box size changes: + +.. code-block:: LAMMPS + + variable x equal "xlo + 0.25*lx" + +Similarly, either of these variable definitions will move the sphere +from an initial position at 2.5 at a constant velocity of 5: + +.. code-block:: LAMMPS + + variable x equal "2.5 + 5*elaplong*dt" + variable x equal vdisplace(2.5,5) + +If a sphere's radius is specified as v_r, then these variable +definitions will grow the size of the sphere at a specified rate. + +.. code-block:: LAMMPS + + variable r0 equal 0.0 + variable rate equal 1.0 + variable r equal "v_r0 + step*dt*v_rate" + +Dump image info +""""""""""""""" + +.. versionadded:: TBD + +Fix graphics is designed to be used with the *fix* keyword of :doc:`dump +image `. The fix will pass geometry information about the +objects listed on the command line to *dump image* so that they are +included in the rendered image. + +The *fflag1* setting of *dump image fix* determines whether cylinder +elements are capped with spheres: 0 means no caps, 1 means the lower +end is capped, 2 means the upper end is capped, and 3 means both ends +are capped. This applies to the *cylinder* object and also to the +body of the *arrow* object and the elements of the *progbar* object. + +The *fflag2* setting allows you to adjust the radius of the rendered +sphere and cylinder items comprising the objects. Since the radius of +these objects is an input parameter for this fix, it is recommended to +set this flag to 0.0. + +.. figure:: JPG/fix-graphics-example.png + :figclass: align-center + + Example of graphics objects rendered with *fix graphics* with *fflag1* setting of 0 (left) and 3 (right) + +These images were created with the following input file: + +.. code-block:: LAMMPS + + units si + region simulation_box block -0.01 0.01 -0.01 0.01 -0.01 0.01 units box + create_box 4 simulation_box + mass * 1 + + variable xpos equal 0.004*sin(PI*step/1000) + variable ypos equal 0.004*cos(PI*step/1000) + variable zpos equal 5.0*v_xpos + variable prog equal (step)/10000.0 + fix gra all graphics 50 sphere 1 v_xpos v_ypos -0.009 0.002 & + sphere 1 0.01 -0.005 0.01 0.005 & + progbar 3 1 z 0.012 -0.012 0.002 0.02 0.0005 v_prog 10 & + cylinder 4 0.01 -0.005 -0.01 0.01 -0.005 0.01 0.005 & + arrow 2 v_xpos v_ypos 0.0 v_xpos v_ypos 0.01 0.0005 0.2 + + dump viz all image 100 myimage2-*.ppm type type size 600 600 zoom 1.24872 & + shiny 0.2 fsaa yes ssao yes 4539 0.6 box no 0.01 axes no 0.5 0.025 & + fix gra type 0 0 view 75 5 + dump_modify viz pad 9 backcolor black acolor 3 gray + + run 2000 + + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +No information about this fix is written to :doc:`binary restart files +`. + +None of the :doc:`fix_modify ` options apply to this fix. + +Restrictions +"""""""""""" + +none + +Related commands +"""""""""""""""" + +none + +Default +""""""" + +none diff --git a/doc/src/fix_indent.rst b/doc/src/fix_indent.rst index 3e269654ac5..1efdc1fb55e 100644 --- a/doc/src/fix_indent.rst +++ b/doc/src/fix_indent.rst @@ -137,9 +137,9 @@ fashion. For the latter, see the *start* and *stop* keywords of the :doc:`thermo_style custom ` for details. For example, if a spherical indenter's x-position is specified as v_x, -then this variable definition will keep it's center at a relative -position in the simulation box, 1/4 of the way from the left edge to -the right edge, even if the box size changes: +then this variable definition will keep its center at a relative +position in the simulation box, 1/4 of the way from the left edge to the +right edge, even if the box size changes: .. code-block:: LAMMPS @@ -201,6 +201,38 @@ contains *xlat*, *ylat*, *zlat* keywords of the :doc:`thermo_style variable k equal 100.0/xlat/xlat fix 1 all indent $k sphere ... +Dump image info +""""""""""""""" + +.. versionadded:: TBD + +Fix indent supports the *fix* keyword of :doc:`dump image `. +The fix will pass geometry information about the indenter to *dump +image* so that the indenter object will be included in the rendered +image. This feature currently only supports spherical, cylindrical, and +planar indenters. Please note, that for :doc:`2d systems `, +a planar indenter rendered as a plane would be invisible and it is thus +rendered as a cylinder. + +The *fflag1* setting of *dump image fix* has no impact on rendering a +spherical indenter or a planar indenter in 3d systems. For a +cylindrical indenter and a planar indenter in 2d systems it determines +whether the cylinder is capped with a sphere at the ends: 0 means no +caps, 1 means the lower end is capped, 2 means the upper end is capped, +and 3 means both ends are capped. + +The *fflag2* setting allows you to adjust the radius of the rendered +object for spherical indenters, cylindrical indenters, and planar +indenters in 2d systems. In many cases you want to use a value < 0 to +reduce the radius of the rendered object so that it does not obscure +atoms close to it. For a planar indenter in 2d systems, it should be +set to a value > 0 or the indenter will not be visible since the +diameter is set internally to zero in that case due to lack of a +suitable heuristic for deriving a meaningful diameter. For a planar +indenter in a 3d system, the *fflag2* value sets the transparency of the +plane. It should be set to a value between 0.0 (invisible) and 1.0 +(fully opaque). + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/doc/src/fix_reaxff_bonds.rst b/doc/src/fix_reaxff_bonds.rst index 748fadf23d4..3cff97ae540 100644 --- a/doc/src/fix_reaxff_bonds.rst +++ b/doc/src/fix_reaxff_bonds.rst @@ -67,6 +67,29 @@ timestep numbers have the same length by adding leading zeroes (e.g. 00010 for a pad value of 5). The default pad value is 0, i.e. no leading zeroes. +.. versionadded:: TBD + +If the filename is "NULL", then no output is created. This can be +useful when using fix *reaxff/bonds* in combination with :doc:`dump +image fix ` keyword to visualize the bonds computed by +the ReaxFF force field. + +Dump image info +""""""""""""""" + +.. versionadded:: TBD + +Fix *reaxff/bonds* supports the *fix* keyword of :doc:`dump image +`. The fix will pass geometry information about the bonds +computed by the :doc:`ReaxFF pair style ` to *dump image* +so that they can be included in the rendered image. + +The *fflag1* setting of *dump image fix* determines whether the bonds +will be capped with spheres (1) or not (0). + +The *fflag2* setting allows to adjust diameter of the cylinders for the +bonds. By default, a diameter of 0.5 length units will be used. + ---------- Restart, fix_modify, output, run start/stop, minimize info diff --git a/doc/src/fix_smd_wall_surface.rst b/doc/src/fix_smd_wall_surface.rst index 203bdb0ca51..53556b940bb 100644 --- a/doc/src/fix_smd_wall_surface.rst +++ b/doc/src/fix_smd_wall_surface.rst @@ -50,6 +50,24 @@ directory. See `this PDF guide `_ to use Smooth Mach Dynamics in LAMMPS. +Dump image info +""""""""""""""" + +.. versionadded:: TBD + +Fix *smd/wall\_surface* supports the *fix* keyword of :doc:`dump image +`. The fix will pass geometry information about the wall +particles to *dump image* so that they be included in the rendered +image. + +The *fflag1* setting of *dump image fix* determines whether the wall will +be rendered as a set of connected triangles (1) or as a mesh of cylinders (2). + +In case of using triangles, the *fflag2* setting determines the +transparency of the triangles and must use a value between 0.0 +(invisible) and 1.0 (fully opaque). If using a mesh of cylinders, the +*fflag2* setting determines the diameter of the cylinders. + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" @@ -59,8 +77,9 @@ minimization. This fix has no outputs. Restrictions """""""""""" -This fix is part of the MACHDYN package. It is only enabled if -LAMMPS was built with that package. See the :doc:`Build package ` page for more info. +This fix is part of the MACHDYN package. It is only enabled if LAMMPS +was built with that package. See the :doc:`Build package +` page for more info. The molecule ID given to the particles created by this fix have to be equal to or larger than 65535. diff --git a/doc/src/fix_wall.rst b/doc/src/fix_wall.rst index c570d85c395..1944f8e5565 100644 --- a/doc/src/fix_wall.rst +++ b/doc/src/fix_wall.rst @@ -497,49 +497,77 @@ if you want the interpolation tables of length Ntable to match exactly what is in the tabulated file (with effectively no preliminary interpolation), you should set Ntable = Nfile. ----------- +----------------- + +Dump image info +""""""""""""""" + +.. versionadded:: TBD + +These wall fixes support the *fix* keyword of :doc:`dump image +`. The fixes will pass geometry information about the walls +to *dump image* so that the walls will be included in the rendered +image. Please note, that for :doc:`2d systems `, a wall +rendered as a plane would be invisible and it is thus rendered as a +cylinder. + +For 2d systems, the *fflag1* setting determines whether the cylinder +representing the wall is capped with a sphere at the ends: 0 means no caps, 1 +means the lower end is capped, 2 means the upper end is capped, and 3 +means both ends are capped. The *fflag2* setting allows to adjust the +radius of the rendered cylinder. It should be set to a value > 0 or the +cylinder will not be visible since the diameter is set internally to +zero due to lack of a suitable heuristic for deriving a meaningful +diameter for all types of walls and unit settings. + +For 3d systems, the *fflag1* setting is ignored, but the *fflag2* +setting determines the transparency of the wall. It must be set to a +value between 0.0 (invisible) and 1.0 (fully opaque). + +----------------- Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" -No information about this fix is written to :doc:`binary restart files -`. +No information about these fixes are written to :doc:`binary restart +files `. -The :doc:`fix_modify ` *energy* option is supported by -this fix to add the energy of interaction between atoms and all the -specified walls to the global potential energy of the system as part -of :doc:`thermodynamic output `. The default setting -for this fix is :doc:`fix_modify energy no `. +The :doc:`fix_modify ` *energy* option is supported by these +fixes to add the energy of interaction between atoms and all the +specified walls to the global potential energy of the system as part of +:doc:`thermodynamic output `. The default setting for +these fixes is :doc:`fix_modify energy no `. The :doc:`fix_modify ` *virial* option is supported by -this fix to add the contribution due to the interaction between atoms +these fixes to add the contribution due to the interaction between atoms and all the specified walls to both the global pressure and per-atom stress of the system via the :doc:`compute pressure ` and :doc:`compute stress/atom ` commands. The former can be accessed by :doc:`thermodynamic output `. The default setting for -this fix is :doc:`fix_modify virial no `. +these fixes is :doc:`fix_modify virial no `. -The :doc:`fix_modify ` *respa* option is supported by this -fix. This allows to set at which level of the :doc:`r-RESPA -` integrator the fix is adding its forces. Default is the +The :doc:`fix_modify ` *respa* option is supported by these +fixes. This allows to set at which level of the :doc:`r-RESPA +` integrator a fix is adding its forces. Default is the outermost level. -This fix computes a global scalar energy and a global vector of forces, -which can be accessed by various :doc:`output commands `. -Note that the scalar energy is the sum of interactions with all defined -walls. If you want the energy on a per-wall basis, you need to use -multiple fix wall commands. The length of the vector is equal to the -number of walls defined by the fix. Each vector value is the normal -force on a specific wall. Note that an outward force on a wall will be -a negative value for *lo* walls and a positive value for *hi* walls. -The scalar and vector values calculated by this fix are "extensive". - -No parameter of this fix can be used with the *start/stop* keywords of -the :doc:`run ` command. - -The forces due to this fix are imposed during an energy minimization, -invoked by the :doc:`minimize ` command. +These fixes compute a global scalar energy and a global vector of +forces, which can be accessed by various :doc:`output commands +`. Note that the scalar energy is the sum of interactions +with all defined walls. If you want the energy on a per-wall basis, you +need to use multiple fix wall commands. The length of the vector is +equal to the number of walls defined by the fix. Each vector value is +the normal force on a specific wall. Note that an outward force on a +wall will be a negative value for *lo* walls and a positive value for +*hi* walls. The scalar and vector values calculated by this fix are +"extensive". + +No parameter of these fixes can be used with the *start/stop* keywords +of the :doc:`run ` command. + +The forces due to these fixes *are* imposed during an energy +minimization, invoked by the :doc:`minimize ` command. .. note:: diff --git a/doc/src/fix_wall_gran.rst b/doc/src/fix_wall_gran.rst index 81a411ffc87..cf3ccd97766 100644 --- a/doc/src/fix_wall_gran.rst +++ b/doc/src/fix_wall_gran.rst @@ -38,15 +38,13 @@ Syntax For *granular*, *fstyle_params* are set using the same syntax as for the *pair_coeff* command of :doc:`pair_style granular ` -* wallstyle = *xplane* or *yplane* or *zplane* or *zcylinder* +* wallstyle = *xplane* or *yplane* or *zplane* * args = list of arguments for a particular style .. parsed-literal:: *xplane* or *yplane* or *zplane* args = lo hi lo,hi = position of lower and upper plane (distance units), either can be NULL) - *zcylinder* args = radius - radius = cylinder radius (distance units) * zero or more keyword/value pairs may be appended to args * keyword = *wiggle* or *shear* or *contacts* or *temperature* @@ -73,7 +71,6 @@ Examples fix 1 all wall/gran hooke 200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0 fix 1 all wall/gran hooke/history 200000.0 NULL 50.0 NULL 0.5 0 zplane 0.0 NULL - fix 2 all wall/gran hooke 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0 fix 3 all wall/gran granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 damping velocity region myBox fix 4 all wall/gran granular jkr 1e5 1500.0 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall region myCone fix 5 all wall/gran granular dmt 1e5 0.2 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall damping tsuji heat 10 region myCone temperature 1.0 @@ -147,29 +144,34 @@ material. system with particles of diameter 1, Kn, Kt, gamma_n, and gamma_s should be set sqrt(2.0) larger than they were previously. -The effective mass *m_eff* in the formulas listed on the :doc:`pair_style granular ` page is the mass of the particle for -particle/wall interactions (mass of wall is infinite). If the +The effective mass *m_eff* in the formulas listed on the +:doc:`pair_style granular ` page is the mass of the particle +for particle/wall interactions (mass of wall is infinite). If the particle is part of a rigid body, its mass is replaced by the mass of -the rigid body in those formulas. This is determined by searching for -a :doc:`fix rigid ` command (or its variants). +the rigid body in those formulas. This is determined by searching for a +:doc:`fix rigid ` command (or its variants). -The *wallstyle* can be planar or cylindrical. The 3 planar options -specify a pair of walls in a dimension. Wall positions are given by -*lo* and *hi*\ . Either of the values can be specified as NULL if a -single wall is desired. For a *zcylinder* wallstyle, the cylinder's -axis is at x = y = 0.0, and the radius of the cylinder is specified. +The *wallstyle* must one of the three planar options. They specify a +pair of walls in a dimension. Wall positions are given by *lo* and +*hi*\ . Either of the values can be specified as NULL if a single wall +is desired. + +.. deprecated:: TBD + +The *zcylinder* wallstyle has been removed. Pleas use :doc:`fix +wall/gran/region ` instead. Optionally, the wall can be moving, if the *wiggle* or *shear* keywords are appended. Both keywords cannot be used together. For the *wiggle* keyword, the wall oscillates sinusoidally, similar to -the oscillations of particles which can be specified by the :doc:`fix move ` command. This is useful in packing simulations of +the oscillations of particles which can be specified by the :doc:`fix +move ` command. This is useful in packing simulations of granular particles. The arguments to the *wiggle* keyword specify a dimension for the motion, as well as it's *amplitude* and *period*\ . Note that if the dimension is in the plane of the wall, this is -effectively a shearing motion. If the dimension is perpendicular to -the wall, it is more of a shaking motion. A *zcylinder* wall can only -be wiggled in the z dimension. +effectively a shearing motion. If the dimension is perpendicular to the +wall, it is more of a shaking motion. Each timestep, the position of a wiggled wall in the appropriate *dim* is set according to this equation: @@ -186,12 +188,11 @@ to the derivative of this expression. For the *shear* keyword, the wall moves continuously in the specified dimension with velocity *vshear*\ . The dimension must be tangential to walls with a planar *wallstyle*, e.g. in the *y* or *z* directions for -an *xplane* wall. For *zcylinder* walls, a dimension of *z* means the -cylinder is moving in the z-direction along it's axis. A dimension of -*x* or *y* means the cylinder is spinning around the z-axis, either in -the clockwise direction for *vshear* > 0 or counter-clockwise for -*vshear* < 0. In this case, *vshear* is the tangential velocity of -the wall at whatever *radius* has been defined. +an *xplane* wall. A dimension of *x* or *y* means the cylinder is +spinning around the z-axis, either in the clockwise direction for +*vshear* > 0 or counter-clockwise for *vshear* < 0. In this case, +*vshear* is the tangential velocity of the wall at whatever *radius* has +been defined. The *temperature* keyword is used to assign a temperature to the wall. The following value can either be a numeric value or an equal-style @@ -261,6 +262,35 @@ No parameter of this fix can be used with the *start/stop* keywords of the :doc:`run ` command. This fix is not invoked during :doc:`energy minimization `. +----------------- + +Dump image info +""""""""""""""" + +.. versionadded:: TBD + +This fix supports the *fix* keyword of :doc:`dump image `. +The fix will pass geometry information about *xplane*\, *yplane*\, and +*zplane* style walls to *dump image* so that the walls will be included +in the rendered image. Please note, that for :doc:`2d systems +`, a wall rendered as a plane would be invisible and it is +thus rendered as a cylinder. + +For 2d systems, the *fflag1* setting determines whether the cylinder +representing the wall is capped with a sphere at the ends: 0 means no caps, 1 +means the lower end is capped, 2 means the upper end is capped, and 3 +means both ends are capped. The *fflag2* setting allows to adjust the +radius of the rendered cylinder. It should be set to a value > 0 or the +cylinder will not be visible since the diameter is set internally to +zero due to lack of a suitable heuristic for deriving a meaningful +diameter for all types of walls and unit settings. + +For 3d systems, the *fflag1* setting is ignored, but the *fflag2* +setting determines the transparency of the wall. It must be set to a +value between 0.0 (invisible) and 1.0 (fully opaque). + +----------------- + Restrictions """""""""""" diff --git a/doc/src/fix_wall_reflect.rst b/doc/src/fix_wall_reflect.rst index c2f4c3a6ab8..01492701772 100644 --- a/doc/src/fix_wall_reflect.rst +++ b/doc/src/fix_wall_reflect.rst @@ -139,6 +139,33 @@ perturbation on the particles: ---------- +Dump image info +""""""""""""""" + +.. versionadded:: TBD + +These wall fixes support the *fix* keyword of :doc:`dump image +`. The fixes will pass geometry information about the walls +to *dump image* so that the walls will be included in the rendered +image. Please note, that for :doc:`2d systems `, a wall +rendered as a plane would be invisible and it is thus rendered as a +cylinder. + +For 2d systems, the *fflag1* setting determines whether the cylinder +representing the wall is capped with a sphere at the ends: 0 means no caps, 1 +means the lower end is capped, 2 means the upper end is capped, and 3 +means both ends are capped. The *fflag2* setting allows to adjust the +radius of the rendered cylinder. It should be set to a value > 0 or the +cylinder will not be visible since the diameter is set internally to +zero due to lack of a suitable heuristic for deriving a meaningful +diameter for all types of walls and unit settings. + +For 3d systems, the *fflag1* setting is ignored, but the *fflag2* +setting determines the transparency of the wall. It must be set to a +value between 0.0 (invisible) and 1.0 (fully opaque). + +---------- + .. include:: accel_styles.rst ---------- diff --git a/doc/src/fix_wall_reflect_stochastic.rst b/doc/src/fix_wall_reflect_stochastic.rst index 5a93950a5da..3f904ca98e3 100644 --- a/doc/src/fix_wall_reflect_stochastic.rst +++ b/doc/src/fix_wall_reflect_stochastic.rst @@ -97,6 +97,30 @@ previously used to define the lattice spacings. ---------- +Dump image info +""""""""""""""" + +.. versionadded:: TBD + +This wall fix supports the *fix* keyword of :doc:`dump image +`. The fix will pass geometry information about the walls +to *dump image* so that the walls will be included in the rendered +image. Please note, that for :doc:`2d systems `, a wall +rendered as a plane would be invisible and it is thus rendered as a +cylinder. + +The *fflag1* setting and the *fflag2* setting of *dump image fix* are +only relevant for 2d systems. The *fflag1* setting determines whether +the cylinder is capped with a sphere at the ends: 0 means no caps, 1 +means the lower end is capped, 2 means the upper end is capped, and 3 +means both ends are capped. The *fflag2* setting allows to adjust the +radius of the rendered cylinder. It should be set to a value > 0 or the +cylinder will not be visible since the diameter is set internally to +zero due to lack of a suitable heuristic for deriving a meaningful +diameter for all types of walls and unit settings. + +---------- + Restrictions """""""""""" diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index d49c8210b7e..f2cc6c1aed2 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -191,6 +191,7 @@ AtomicPairStyle atomID atomistic atomtypes +atrans attogram attograms attrac @@ -217,6 +218,7 @@ awpmd AWPMD ax Axel +axestrans Axilrod Ay ay @@ -384,6 +386,7 @@ Bourne boxcolor boxhi boxlo +boxtrans boxxhi boxxlo boxyhi @@ -415,6 +418,7 @@ Bruskin Brusselle Bryantsev Btarget +btrans btype buckPlusAttr buf @@ -1736,6 +1740,7 @@ isokinetic isomorphism isothermal isotropically +isovalue isovolume Isralewitz iter @@ -3758,6 +3763,7 @@ strstr Stukowski Su subbox +subboxtrans Subclassed subcutoff subcycle diff --git a/examples/indent/in.indent b/examples/indent/in.indent index ba19e93cf51..6338d86fa5b 100644 --- a/examples/indent/in.indent +++ b/examples/indent/in.indent @@ -1,67 +1,71 @@ # 2d indenter simulation -dimension 2 -boundary p s p +dimension 2 +boundary p s p -atom_style atomic -neighbor 0.3 bin -neigh_modify delay 5 +atom_style atomic +neighbor 0.3 bin +neigh_modify delay 5 # create geometry -lattice hex 0.9 -region box block 0 20 0 10 -0.25 0.25 -create_box 2 box -create_atoms 1 box +lattice hex 0.9 +region box block 0 20 0 10 -0.25 0.25 +create_box 2 box +create_atoms 1 box -mass 1 1.0 -mass 2 1.0 +mass 1 1.0 +mass 2 1.0 # LJ potentials -pair_style lj/cut 2.5 -pair_coeff * * 1.0 1.0 2.5 +pair_style lj/cut 2.5 +pair_coeff * * 1.0 1.0 2.5 # define groups -region 1 block INF INF INF 1.25 INF INF -group lower region 1 -group mobile subtract all lower -set group lower type 2 +region 1 block INF INF INF 1.25 INF INF +group lower region 1 +group mobile subtract all lower +set group lower type 2 # initial velocities -compute new mobile temp -velocity mobile create 0.2 482748 temp new -fix 1 all nve -fix 2 lower setforce 0.0 0.0 0.0 -fix 3 all temp/rescale 100 0.1 0.1 0.01 1.0 +compute new mobile temp +velocity mobile create 0.2 482748 temp new +fix 1 all nve +fix 2 lower setforce 0.0 0.0 0.0 +fix 3 all temp/rescale 100 0.1 0.1 0.01 1.0 # run with indenter -timestep 0.003 -variable k equal 1000.0/xlat +timestep 0.003 +variable k equal 1000.0/xlat variable y equal "13.0*ylat - step*dt*0.02*ylat" -fix 4 all indent $k sphere 10 v_y 0 5.0 -fix 5 all enforce2d +fix 4 all indent $k sphere 10 v_y 0 5.0 +fix 5 all enforce2d -thermo 1000 -thermo_modify temp new +thermo 1000 +thermo_modify temp new -#dump 1 all atom 250 dump.indent +#dump 1 all atom 250 dump.indent -#dump 2 all image 1000 image.*.jpg type type & -# zoom 1.6 adiam 1.5 -#dump_modify 2 pad 5 +#dump 2 all image 1000 image.*.jpg type type & +# zoom 1.6 adiam 1.5 fix 4 type 0.0 -2.0 +#dump_modify 2 pad 5 fcolor purple -#dump 3 all movie 1000 movie.mpg type type & -# zoom 1.6 adiam 1.5 -#dump_modify 3 pad 5 +#dump 3 all movie 1000 movie.mpg type type & +# zoom 1.6 adiam 1.5 +#dump_modify 3 pad 5 -run 30000 +run 30000 # run without indenter unfix 4 -run 30000 +#undump 2 +#dump 2 all image 1000 image.*.jpg type type & +# zoom 1.6 adiam 1.5 +#dump_modify 2 pad 5 +run 30000 diff --git a/src/BODY/body_nparticle.cpp b/src/BODY/body_nparticle.cpp index 2d4c06ae9e5..0e1a8fc30dc 100644 --- a/src/BODY/body_nparticle.cpp +++ b/src/BODY/body_nparticle.cpp @@ -16,6 +16,7 @@ #include "atom.h" #include "atom_vec_body.h" +#include "dump_image.h" #include "error.h" #include "math_extra.h" #include "math_eigen.h" @@ -27,7 +28,6 @@ using namespace LAMMPS_NS; static constexpr double EPSILON = 1.0e-7; -enum{SPHERE,LINE,TRI}; // also in DumpImage /* ---------------------------------------------------------------------- */ @@ -355,7 +355,7 @@ int BodyNparticle::image(int ibonus, double flag1, double /*flag2*/, int n = bonus->ivalue[0]; for (int i = 0; i < n; i++) { - imflag[i] = SPHERE; + imflag[i] = DumpImage::SPHERE; MathExtra::quat_to_mat(bonus->quat,p); MathExtra::matvec(p,&bonus->dvalue[3*i],imdata[i]); diff --git a/src/BODY/body_rounded_polygon.cpp b/src/BODY/body_rounded_polygon.cpp index 366db6264f1..3fc366440fd 100644 --- a/src/BODY/body_rounded_polygon.cpp +++ b/src/BODY/body_rounded_polygon.cpp @@ -21,6 +21,7 @@ #include "atom.h" #include "atom_vec_body.h" #include "domain.h" +#include "dump_image.h" #include "error.h" #include "math_extra.h" #include "math_eigen.h" @@ -33,7 +34,6 @@ using namespace LAMMPS_NS; static constexpr double EPSILON = 1.0e-7; -enum{SPHERE,LINE}; // also in DumpImage /* ---------------------------------------------------------------------- */ @@ -510,7 +510,7 @@ int BodyRoundedPolygon::image(int ibonus, double flag1, double /*flag2*/, if (n == 1) { for (int i = 0; i < n; i++) { - imflag[i] = SPHERE; + imflag[i] = DumpImage::SPHERE; MathExtra::quat_to_mat(bonus->quat,p); MathExtra::matvec(p,&bonus->dvalue[3*i],imdata[i]); @@ -528,7 +528,7 @@ int BodyRoundedPolygon::image(int ibonus, double flag1, double /*flag2*/, // first end pt of each line for (int i = 0; i < n; i++) { - imflag[i] = LINE; + imflag[i] = DumpImage::LINE; MathExtra::quat_to_mat(bonus->quat,p); MathExtra::matvec(p,&bonus->dvalue[3*i],imdata[i]); diff --git a/src/BODY/body_rounded_polyhedron.cpp b/src/BODY/body_rounded_polyhedron.cpp index bd16dac96c3..30b757eb2bc 100644 --- a/src/BODY/body_rounded_polyhedron.cpp +++ b/src/BODY/body_rounded_polyhedron.cpp @@ -20,6 +20,7 @@ #include "atom.h" #include "atom_vec_body.h" +#include "dump_image.h" #include "error.h" #include "math_extra.h" #include "math_eigen.h" @@ -34,8 +35,6 @@ using namespace LAMMPS_NS; static constexpr double EPSILON = 1.0e-7; static constexpr int MAX_FACE_SIZE = 4; // maximum number of vertices per face (for now) -enum { SPHERE, LINE }; // also in DumpImage - /* ---------------------------------------------------------------------- */ BodyRoundedPolyhedron::BodyRoundedPolyhedron(LAMMPS *lmp, int narg, char **arg) : @@ -615,7 +614,7 @@ int BodyRoundedPolyhedron::image(int ibonus, double flag1, double /*flag2*/, if (nvertices == 1) { // spheres for (int i = 0; i < nvertices; i++) { - imflag[i] = SPHERE; + imflag[i] = DumpImage::SPHERE; MathExtra::quat_to_mat(bonus->quat,p); MathExtra::matvec(p,&bonus->dvalue[3*i],imdata[i]); @@ -637,7 +636,7 @@ int BodyRoundedPolyhedron::image(int ibonus, double flag1, double /*flag2*/, int pt1, pt2; for (int i = 0; i < nedges; i++) { - imflag[i] = LINE; + imflag[i] = DumpImage::LINE; pt1 = static_cast(edge_ends[2*i]); pt2 = static_cast(edge_ends[2*i+1]); diff --git a/src/EXTRA-FIX/fix_wall_reflect_stochastic.cpp b/src/EXTRA-FIX/fix_wall_reflect_stochastic.cpp index 42bb1980495..d07a25409a4 100644 --- a/src/EXTRA-FIX/fix_wall_reflect_stochastic.cpp +++ b/src/EXTRA-FIX/fix_wall_reflect_stochastic.cpp @@ -23,6 +23,7 @@ #include "comm.h" #include "domain.h" #include "error.h" +#include "fix_wall.h" #include "force.h" #include "lattice.h" #include "random_mars.h" @@ -79,12 +80,12 @@ FixWallReflectStochastic(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); int newwall; - if (strcmp(arg[iarg],"xlo") == 0) newwall = XLO; - else if (strcmp(arg[iarg],"xhi") == 0) newwall = XHI; - else if (strcmp(arg[iarg],"ylo") == 0) newwall = YLO; - else if (strcmp(arg[iarg],"yhi") == 0) newwall = YHI; - else if (strcmp(arg[iarg],"zlo") == 0) newwall = ZLO; - else if (strcmp(arg[iarg],"zhi") == 0) newwall = ZHI; + if (strcmp(arg[iarg],"xlo") == 0) newwall = FixWall::XLO; + else if (strcmp(arg[iarg],"xhi") == 0) newwall = FixWall::XHI; + else if (strcmp(arg[iarg],"ylo") == 0) newwall = FixWall::YLO; + else if (strcmp(arg[iarg],"yhi") == 0) newwall = FixWall::YHI; + else if (strcmp(arg[iarg],"zlo") == 0) newwall = FixWall::ZLO; + else if (strcmp(arg[iarg],"zhi") == 0) newwall = FixWall::ZHI; for (int m = 0; (m < nwall) && (m < 6); m++) if (newwall == wallwhich[m]) @@ -139,19 +140,19 @@ FixWallReflectStochastic(LAMMPS *lmp, int narg, char **arg) : if (nwall == 0) error->all(FLERR,"Illegal fix wall command"); for (int m = 0; m < nwall; m++) { - if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic) + if ((wallwhich[m] == FixWall::XLO || wallwhich[m] == FixWall::XHI) && domain->xperiodic) error->all(FLERR,"Cannot use fix wall/reflect/stochastic " "in periodic dimension"); - if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic) + if ((wallwhich[m] == FixWall::YLO || wallwhich[m] == FixWall::YHI) && domain->yperiodic) error->all(FLERR,"Cannot use fix wall/reflect/stochastic " "in periodic dimension"); - if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic) + if ((wallwhich[m] == FixWall::ZLO || wallwhich[m] == FixWall::ZHI) && domain->zperiodic) error->all(FLERR,"Cannot use fix wall/reflect/stochastic " "in periodic dimension"); } for (int m = 0; m < nwall; m++) - if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2) + if ((wallwhich[m] == FixWall::ZLO || wallwhich[m] == FixWall::ZHI) && domain->dimension == 2) error->all(FLERR, "Cannot use fix wall/reflect/stochastic zlo/zhi " "for a 2d simulation"); @@ -171,8 +172,8 @@ FixWallReflectStochastic(LAMMPS *lmp, int narg, char **arg) : for (int m = 0; m < nwall; m++) { if (wallstyle[m] != CONSTANT) continue; - if (wallwhich[m] < YLO) coord0[m] *= xscale; - else if (wallwhich[m] < ZLO) coord0[m] *= yscale; + if (wallwhich[m] < FixWall::YLO) coord0[m] *= xscale; + else if (wallwhich[m] < FixWall::ZLO) coord0[m] *= yscale; else coord0[m] *= zscale; } } diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index d58ce011f96..5dd8d6a9fc1 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -23,7 +23,9 @@ #include "granular_model.h" #include "gran_sub_mod.h" #include "domain.h" +#include "dump_image.h" #include "error.h" +#include "fix_wall.h" #include "input.h" #include "math_const.h" #include "math_extra.h" @@ -47,14 +49,14 @@ static constexpr double BIG = 1.0e20; // XYZ PLANE need to be 0,1,2 -enum {NOSTYLE=-1,XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER,REGION}; +enum {NOSTYLE=-1,XPLANE=0,YPLANE=1,ZPLANE=2,REGION}; enum {NONE,CONSTANT,EQUAL}; /* ---------------------------------------------------------------------- */ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), idregion(nullptr), tstr(nullptr), history_one(nullptr), - fix_rigid(nullptr), mass_rigid(nullptr) + fix_rigid(nullptr), mass_rigid(nullptr), imgobjs(nullptr), imgparms(nullptr) { if (narg < 4) utils::missing_cmd_args(FLERR,"fix wall/gran", error); @@ -128,41 +130,68 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : // wallstyle args + numwalls = 0; if (iarg >= narg) error->all(FLERR, "Illegal fix wall/gran command"); if (strcmp(arg[iarg],"xplane") == 0) { if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/gran command"); wallstyle = XPLANE; - if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG; - else lo = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG; - else hi = utils::numeric(FLERR,arg[iarg+2],false,lmp); + numwalls = 0; + if (strcmp(arg[iarg+1],"NULL") == 0) { + lo = -BIG; + } else { + lo = utils::numeric(FLERR,arg[iarg+1],false,lmp); + ++numwalls; + } + if (strcmp(arg[iarg+2],"NULL") == 0) { + hi = BIG; + } else { + hi = utils::numeric(FLERR,arg[iarg+2],false,lmp); + ++numwalls; + } iarg += 3; } else if (strcmp(arg[iarg],"yplane") == 0) { if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/gran command"); wallstyle = YPLANE; - if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG; - else lo = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG; - else hi = utils::numeric(FLERR,arg[iarg+2],false,lmp); + numwalls = 0; + if (strcmp(arg[iarg+1],"NULL") == 0) { + lo = -BIG; + } else { + lo = utils::numeric(FLERR,arg[iarg+1],false,lmp); + ++numwalls; + } + if (strcmp(arg[iarg+2],"NULL") == 0) { + hi = BIG; + } else { + hi = utils::numeric(FLERR,arg[iarg+2],false,lmp); + ++numwalls; + } iarg += 3; } else if (strcmp(arg[iarg],"zplane") == 0) { if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/gran command"); wallstyle = ZPLANE; - if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG; - else lo = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG; - else hi = utils::numeric(FLERR,arg[iarg+2],false,lmp); + numwalls = 0; + if (strcmp(arg[iarg+1],"NULL") == 0) { + lo = -BIG; + } else { + lo = utils::numeric(FLERR,arg[iarg+1],false,lmp); + ++numwalls; + } + if (strcmp(arg[iarg+2],"NULL") == 0) { + hi = BIG; + } else { + hi = utils::numeric(FLERR,arg[iarg+2],false,lmp); + ++numwalls; + } iarg += 3; } else if (strcmp(arg[iarg],"zcylinder") == 0) { - if (narg < iarg+2) error->all(FLERR,"Illegal fix wall/gran command"); - wallstyle = ZCYLINDER; - lo = hi = 0.0; - cylradius = utils::numeric(FLERR,arg[iarg+1],false,lmp); iarg += 2; + error->all(FLERR, iarg, "The zcylinder keyword has been removed. " + "Please use fix wall/gran/region instead."); } else if (strcmp(arg[iarg],"region") == 0) { if (narg < iarg+2) error->all(FLERR,"Illegal fix wall/gran command"); wallstyle = REGION; + numwalls = 0; delete[] idregion; idregion = utils::strdup(arg[iarg+1]); iarg += 2; @@ -227,13 +256,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Cannot use wall in periodic dimension"); if (wallstyle == ZPLANE && domain->zperiodic) error->all(FLERR,"Cannot use wall in periodic dimension"); - if (wallstyle == ZCYLINDER && (domain->xperiodic || domain->yperiodic)) - error->all(FLERR,"Cannot use wall in periodic dimension"); if (wiggle && wshear) error->all(FLERR,"Cannot wiggle and shear fix wall/gran"); - if (wiggle && wallstyle == ZCYLINDER && axis != 2) - error->all(FLERR,"Invalid wiggle direction for fix wall/gran"); if (wshear && wallstyle == XPLANE && axis == 0) error->all(FLERR,"Invalid shear direction for fix wall/gran"); if (wshear && wallstyle == YPLANE && axis == 1) @@ -271,6 +296,29 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : } time_origin = update->ntimestep; + + // for rendering walls with dump image. + if ((wallstyle == XPLANE) || (wallstyle == YPLANE) || (wallstyle == ZPLANE)) { + if (domain->dimension == 2) { + // one cylinder object per wall to draw in 2d + memory->create(imgobjs, numwalls, "fix_wall:imgobjs"); + memory->create(imgparms, numwalls, 8, "fix_wall:imgparms"); + for (int m = 0; m < numwalls; ++m) { + imgobjs[m] = DumpImage::CYLINDER; + imgparms[m][0] = 1; // use color of first atom type by default + } + } else { + // two triangle objects per wall to draw in 3d + memory->create(imgobjs, 2 * numwalls, "fix_wall:imgobjs"); + memory->create(imgparms, 2 * numwalls, 10, "fix_wall:imgparms"); + for (int m = 0; m < numwalls; ++m) { + imgobjs[2 * m] = DumpImage::TRIANGLE; + imgobjs[2 * m + 1] = DumpImage::TRIANGLE; + imgparms[2 * m][0] = 1; // use color of first atom type by default + imgparms[2 * m + 1][0] = 1; // use color of first atom type by default + } + } + } } /* ---------------------------------------------------------------------- */ @@ -291,6 +339,9 @@ FixWallGran::~FixWallGran() delete[] idregion; memory->destroy(history_one); memory->destroy(mass_rigid); + + memory->destroy(imgobjs); + memory->destroy(imgparms); } /* ---------------------------------------------------------------------- */ @@ -475,23 +526,6 @@ void FixWallGran::post_force(int /*vflag*/) del2 = whi - x[i][2]; if (del1 < del2) dz = del1; else dz = -del2; - } else if (wallstyle == ZCYLINDER) { - delxy = sqrt(x[i][0] * x[i][0] + x[i][1] * x[i][1]); - delr = cylradius - delxy; - if (delr > radius[i]) { - dz = cylradius; - rwall = 0.0; - } else { - dx = -delr / delxy * x[i][0]; - dy = -delr / delxy * x[i][1]; - // rwall = -2r_c if inside cylinder, 2r_c outside - rwall = (delxy < cylradius) ? -2 * cylradius : 2 * cylradius; - if (wshear && axis != 2) { - vwall[0] += vshear * x[i][1] / delxy; - vwall[1] += -vshear * x[i][0] / delxy; - vwall[2] = 0.0; - } - } } // Reset model and copy initial geometric data @@ -553,6 +587,37 @@ void FixWallGran::post_force(int /*vflag*/) array_atom[i][8 + n] = model->svector[n]; } } + + // update wall image information + int m = 0; + if (wallstyle == XPLANE) { + if (lo != -BIG) { + FixWall::update_image_plane(m, FixWall::XLO, wlo, imgparms, domain); + ++m; + } + if (hi != BIG) { + FixWall::update_image_plane(m, FixWall::XHI, whi, imgparms, domain); + ++m; + } + } else if (wallstyle == YPLANE) { + if (lo != -BIG) { + FixWall::update_image_plane(m, FixWall::YLO, wlo, imgparms, domain); + ++m; + } + if (hi != BIG) { + FixWall::update_image_plane(m, FixWall::YHI, whi, imgparms, domain); + ++m; + } + } else if (wallstyle == ZPLANE) { + if (lo != -BIG) { + FixWall::update_image_plane(m, FixWall::ZLO, wlo, imgparms, domain); + ++m; + } + if (hi != BIG) { + FixWall::update_image_plane(m, FixWall::ZHI, whi, imgparms, domain); + ++m; + } + } } /* ---------------------------------------------------------------------- */ @@ -731,3 +796,29 @@ void FixWallGran::reset_dt() dt = update->dt; model->dt = dt; } + +/* ---------------------------------------------------------------------- + provide graphics information to dump image to render wall as plane + data has been copied to dedicated storage during fix indent execution +------------------------------------------------------------------------- */ + +int FixWallGran::image(int *&objs, double **&parms) +{ + objs = imgobjs; + parms = imgparms; + switch (wallstyle) { + case XPLANE: // fallthrough + case YPLANE: // fallthrough + case ZPLANE: + if (domain->dimension == 2) { + return numwalls; + } else { + return 2 * numwalls; + } + break; + case REGION: // can visualize region directly + return 0; + break; + } + return 0; +} diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index e3d35fcc36a..abb769ee5bb 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -50,6 +50,8 @@ class FixWallGran : public Fix { int maxsize_restart() override; void reset_dt() override; + int image(int *&, double **&) override; + // for granular model choices class Granular_NS::GranularModel *model; @@ -58,7 +60,7 @@ class FixWallGran : public Fix { int nlevels_respa; bigint time_origin; - double lo, hi, cylradius; + double lo, hi; double amplitude, period, omega, vshear; double dt; double Twall; @@ -82,6 +84,12 @@ class FixWallGran : public Fix { double *mass_rigid; // rigid mass for owned+ghost atoms int nmax; // allocated size of mass_rigid + // dump image data + + int numwalls; + int *imgobjs; + double **imgparms; + // store particle interactions int nsvector; diff --git a/src/KOKKOS/fix_wall_reflect_kokkos.cpp b/src/KOKKOS/fix_wall_reflect_kokkos.cpp index 21a7b84aef1..4e5262fcb58 100644 --- a/src/KOKKOS/fix_wall_reflect_kokkos.cpp +++ b/src/KOKKOS/fix_wall_reflect_kokkos.cpp @@ -17,6 +17,7 @@ #include "atom_kokkos.h" #include "atom_masks.h" #include "input.h" +#include "fix_wall.h" #include "modify.h" #include "update.h" #include "variable.h" @@ -58,14 +59,17 @@ void FixWallReflectKokkos::post_integrate() for (int m = 0; m < nwall; m++) { if (wallstyle[m] == VARIABLE) { coord = input->variable->compute_equal(varindex[m]); - if (wallwhich[m] < YLO) coord *= xscale; - else if (wallwhich[m] < ZLO) coord *= yscale; + if (wallwhich[m] < FixWall::YLO) coord *= xscale; + else if (wallwhich[m] < FixWall::ZLO) coord *= yscale; else coord *= zscale; } else coord = coord0[m]; dim = wallwhich[m] / 2; side = wallwhich[m] % 2; + // record wall graphics objects for dump image + FixWall::update_image_plane(m, wallwhich[m], coord, imgparms, domain); + copymode = 1; Kokkos::parallel_for(Kokkos::RangePolicy(0,nlocal),*this); copymode = 0; diff --git a/src/MACHDYN/fix_smd_wall_surface.cpp b/src/MACHDYN/fix_smd_wall_surface.cpp index 25e76e1dab2..8510f722cb2 100644 --- a/src/MACHDYN/fix_smd_wall_surface.cpp +++ b/src/MACHDYN/fix_smd_wall_surface.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -22,12 +21,14 @@ #include "atom_vec.h" #include "comm.h" #include "domain.h" +#include "dump_image.h" #include "error.h" +#include "memory.h" #include "text_file_reader.h" +#include #include #include -#include using namespace LAMMPS_NS; using namespace FixConst; @@ -39,65 +40,63 @@ static constexpr double EPSILON = 1.0e-6; /* ---------------------------------------------------------------------- */ FixSMDWallSurface::FixSMDWallSurface(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) { - - restart_global = 0; - restart_peratom = 0; - first = 1; - - //atom->add_callback(Atom::GROW); - //atom->add_callback(Atom::RESTART); - - if (narg != 6) - error->all(FLERR, "Illegal number of arguments for fix smd/wall_surface"); - - filename = strdup(arg[3]); - wall_particle_type = utils::inumeric(FLERR,arg[4],false,lmp); - wall_molecule_id = utils::inumeric(FLERR,arg[5],false,lmp); - if (wall_molecule_id < 65535) { - error->one(FLERR, "wall molcule id must be >= 65535\n"); - } - - if (comm->me == 0) { - printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n"); - printf("fix smd/wall_surface reads trianglulated surface from file: %s\n", filename); - printf("fix smd/wall_surface has particle type %d \n", wall_particle_type); - printf("fix smd/wall_surface has molecule id %d \n", wall_molecule_id); - printf(">>========>>========>>========>>========>>========>>========>>========>>========\n"); - } + Fix(lmp, narg, arg), imgobjs(nullptr), imgparms(nullptr) +{ + + restart_global = 0; + restart_peratom = 0; + first = 1; + + if (narg != 6) error->all(FLERR, "Illegal number of arguments for fix smd/wall_surface"); + + filename = strdup(arg[3]); + wall_particle_type = utils::inumeric(FLERR, arg[4], false, lmp); + wall_molecule_id = utils::inumeric(FLERR, arg[5], false, lmp); + if (wall_molecule_id < 65535) { error->one(FLERR, "wall molcule id must be >= 65535\n"); } + + if (comm->me == 0) { + printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n"); + printf("fix smd/wall_surface reads trianglulated surface from file: %s\n", filename); + printf("fix smd/wall_surface has particle type %d \n", wall_particle_type); + printf("fix smd/wall_surface has molecule id %d \n", wall_molecule_id); + printf(">>========>>========>>========>>========>>========>>========>>========>>========\n"); + } } /* ---------------------------------------------------------------------- */ -FixSMDWallSurface::~FixSMDWallSurface() { - free(filename); - filename = nullptr; - // unregister this fix so atom class doesn't invoke it any more +FixSMDWallSurface::~FixSMDWallSurface() +{ + free(filename); + filename = nullptr; - //atom->delete_callback(id,Atom::GROW); - //atom->delete_callback(id,Atom::RESTART); + // clean up dump image data + memory->destroy(imgobjs); + memory->destroy(imgparms); } /* ---------------------------------------------------------------------- */ -int FixSMDWallSurface::setmask() { - int mask = 0; - return mask; +int FixSMDWallSurface::setmask() +{ + int mask = 0; + return mask; } /* ---------------------------------------------------------------------- */ -void FixSMDWallSurface::init() { - if (!first) - return; +void FixSMDWallSurface::init() +{ + if (!first) return; } /* ---------------------------------------------------------------------- For minimization: setup as with dynamics ------------------------------------------------------------------------- */ -void FixSMDWallSurface::min_setup(int vflag) { - setup(vflag); +void FixSMDWallSurface::min_setup(int vflag) +{ + setup(vflag); } /* ---------------------------------------------------------------------- @@ -105,92 +104,81 @@ void FixSMDWallSurface::min_setup(int vflag) { must be done in setup (not init) since fix init comes before neigh init ------------------------------------------------------------------------- */ -void FixSMDWallSurface::setup(int /*vflag*/) { - - if (!first) - return; - first = 0; - - // set bounds for my proc - // if periodic and I am lo/hi proc, adjust bounds by EPSILON - // ensures all data atoms will be owned even with round-off - - int triclinic = domain->triclinic; - - double epsilon[3]; - if (triclinic) - epsilon[0] = epsilon[1] = epsilon[2] = EPSILON; - else { - epsilon[0] = domain->prd[0] * EPSILON; - epsilon[1] = domain->prd[1] * EPSILON; - epsilon[2] = domain->prd[2] * EPSILON; - } - - if (triclinic == 0) { - sublo[0] = domain->sublo[0]; - subhi[0] = domain->subhi[0]; - sublo[1] = domain->sublo[1]; - subhi[1] = domain->subhi[1]; - sublo[2] = domain->sublo[2]; - subhi[2] = domain->subhi[2]; - } else { - sublo[0] = domain->sublo_lamda[0]; - subhi[0] = domain->subhi_lamda[0]; - sublo[1] = domain->sublo_lamda[1]; - subhi[1] = domain->subhi_lamda[1]; - sublo[2] = domain->sublo_lamda[2]; - subhi[2] = domain->subhi_lamda[2]; - } - - if (comm->layout != Comm::LAYOUT_TILED) { - if (domain->xperiodic) { - if (comm->myloc[0] == 0) - sublo[0] -= epsilon[0]; - if (comm->myloc[0] == comm->procgrid[0] - 1) - subhi[0] += epsilon[0]; - } - if (domain->yperiodic) { - if (comm->myloc[1] == 0) - sublo[1] -= epsilon[1]; - if (comm->myloc[1] == comm->procgrid[1] - 1) - subhi[1] += epsilon[1]; - } - if (domain->zperiodic) { - if (comm->myloc[2] == 0) - sublo[2] -= epsilon[2]; - if (comm->myloc[2] == comm->procgrid[2] - 1) - subhi[2] += epsilon[2]; - } - - } else { - if (domain->xperiodic) { - if (comm->mysplit[0][0] == 0.0) - sublo[0] -= epsilon[0]; - if (comm->mysplit[0][1] == 1.0) - subhi[0] += epsilon[0]; - } - if (domain->yperiodic) { - if (comm->mysplit[1][0] == 0.0) - sublo[1] -= epsilon[1]; - if (comm->mysplit[1][1] == 1.0) - subhi[1] += epsilon[1]; - } - if (domain->zperiodic) { - if (comm->mysplit[2][0] == 0.0) - sublo[2] -= epsilon[2]; - if (comm->mysplit[2][1] == 1.0) - subhi[2] += epsilon[2]; - } - } - - read_triangles(0); +void FixSMDWallSurface::setup(int /*vflag*/) +{ + + if (!first) return; + first = 0; + + // set bounds for my proc + // if periodic and I am lo/hi proc, adjust bounds by EPSILON + // ensures all data atoms will be owned even with round-off + + int triclinic = domain->triclinic; + + double epsilon[3]; + if (triclinic) + epsilon[0] = epsilon[1] = epsilon[2] = EPSILON; + else { + epsilon[0] = domain->prd[0] * EPSILON; + epsilon[1] = domain->prd[1] * EPSILON; + epsilon[2] = domain->prd[2] * EPSILON; + } + + if (triclinic == 0) { + sublo[0] = domain->sublo[0]; + subhi[0] = domain->subhi[0]; + sublo[1] = domain->sublo[1]; + subhi[1] = domain->subhi[1]; + sublo[2] = domain->sublo[2]; + subhi[2] = domain->subhi[2]; + } else { + sublo[0] = domain->sublo_lamda[0]; + subhi[0] = domain->subhi_lamda[0]; + sublo[1] = domain->sublo_lamda[1]; + subhi[1] = domain->subhi_lamda[1]; + sublo[2] = domain->sublo_lamda[2]; + subhi[2] = domain->subhi_lamda[2]; + } + + if (comm->layout != Comm::LAYOUT_TILED) { + if (domain->xperiodic) { + if (comm->myloc[0] == 0) sublo[0] -= epsilon[0]; + if (comm->myloc[0] == comm->procgrid[0] - 1) subhi[0] += epsilon[0]; + } + if (domain->yperiodic) { + if (comm->myloc[1] == 0) sublo[1] -= epsilon[1]; + if (comm->myloc[1] == comm->procgrid[1] - 1) subhi[1] += epsilon[1]; + } + if (domain->zperiodic) { + if (comm->myloc[2] == 0) sublo[2] -= epsilon[2]; + if (comm->myloc[2] == comm->procgrid[2] - 1) subhi[2] += epsilon[2]; + } + + } else { + if (domain->xperiodic) { + if (comm->mysplit[0][0] == 0.0) sublo[0] -= epsilon[0]; + if (comm->mysplit[0][1] == 1.0) subhi[0] += epsilon[0]; + } + if (domain->yperiodic) { + if (comm->mysplit[1][0] == 0.0) sublo[1] -= epsilon[1]; + if (comm->mysplit[1][1] == 1.0) subhi[1] += epsilon[1]; + } + if (domain->zperiodic) { + if (comm->mysplit[2][0] == 0.0) sublo[2] -= epsilon[2]; + if (comm->mysplit[2][1] == 1.0) subhi[2] += epsilon[2]; + } + } + + read_triangles(0); } /* ---------------------------------------------------------------------- size of atom nlocal's restart data ------------------------------------------------------------------------- */ -void FixSMDWallSurface::read_triangles(int pass) { +void FixSMDWallSurface::read_triangles(int pass) +{ double coord[3]; int nlocal_previous = atom->nlocal; @@ -201,8 +189,7 @@ void FixSMDWallSurface::read_triangles(int pass) { Vector3d normal, center; FILE *fp = fopen(filename, "r"); - if (fp == nullptr) - error->one(FLERR, "Cannot open file {}: {}", filename, utils::getsyserror()); + if (fp == nullptr) error->one(FLERR, "Cannot open file {}: {}", filename, utils::getsyserror()); if (comm->me == 0) { utils::logmesg(lmp, "\n>>========>>========>>========>>========>>========>>========\n"); @@ -216,43 +203,43 @@ void FixSMDWallSurface::read_triangles(int pass) { try { char *line = reader.next_line(); if (!line || !utils::strmatch(line, "^solid")) - throw TokenizerException("Invalid triangles file format",""); + throw TokenizerException("Invalid triangles file format", ""); if (comm->me == 0) - utils::logmesg(lmp, " reading STL object '{}' from {}\n", utils::trim(line+6), filename); + utils::logmesg(lmp, " reading STL object '{}' from {}\n", utils::trim(line + 6), filename); - while((line = reader.next_line())) { + while ((line = reader.next_line())) { // next line is facet line with 5 words auto values = utils::split_words(line); // otherwise stop reading - if ((values.size() != 5) || !utils::strmatch(values[0],"^facet")) break; + if ((values.size() != 5) || !utils::strmatch(values[0], "^facet")) break; normal << utils::numeric(FLERR, values[2], false, lmp), - utils::numeric(FLERR, values[3], false, lmp), - utils::numeric(FLERR, values[4], false, lmp); + utils::numeric(FLERR, values[3], false, lmp), + utils::numeric(FLERR, values[4], false, lmp); line = reader.next_line(2); if (!line || !utils::strmatch(line, "^ *outer *loop")) - throw TokenizerException("Error reading outer loop",""); + throw TokenizerException("Error reading outer loop", ""); for (int k = 0; k < 3; ++k) { line = reader.next_line(4); values = utils::split_words(line); - if ((values.size() != 4) || !utils::strmatch(values[0],"^vertex")) - throw TokenizerException("Error reading vertex",""); + if ((values.size() != 4) || !utils::strmatch(values[0], "^vertex")) + throw TokenizerException("Error reading vertex", ""); vert[k] << utils::numeric(FLERR, values[1], false, lmp), - utils::numeric(FLERR, values[2], false, lmp), - utils::numeric(FLERR, values[3], false, lmp); + utils::numeric(FLERR, values[2], false, lmp), + utils::numeric(FLERR, values[3], false, lmp); } line = reader.next_line(1); if (!line || !utils::strmatch(line, "^ *endloop")) - throw TokenizerException("Error reading endloop",""); + throw TokenizerException("Error reading endloop", ""); line = reader.next_line(1); if (!line || !utils::strmatch(line, "^ *endfacet")) - throw TokenizerException("Error reading endfacet",""); + throw TokenizerException("Error reading endfacet", ""); // now we have a normal and three vertices ... proceed with adding triangle @@ -271,8 +258,8 @@ void FixSMDWallSurface::read_triangles(int pass) { * ... radius is the mmaximal distance from triangle center to all vertices */ - if (center(0) >= sublo[0] && center(0) < subhi[0] && center(1) >= sublo[1] && center(1) < subhi[1] && center(2) >= sublo[2] - && center(2) < subhi[2]) { + if (center(0) >= sublo[0] && center(0) < subhi[0] && center(1) >= sublo[1] && + center(1) < subhi[1] && center(2) >= sublo[2] && center(2) < subhi[2]) { coord[0] = center(0); coord[1] = center(1); @@ -291,8 +278,8 @@ void FixSMDWallSurface::read_triangles(int pass) { double **smd_data_9 = atom->smd_data_9; double **x0 = atom->x0; - radius[ilocal] = r; //ilocal; - contact_radius[ilocal] = r; //ilocal; + radius[ilocal] = r; //ilocal; + contact_radius[ilocal] = r; //ilocal; mol[ilocal] = wall_molecule_id; type[ilocal] = wall_particle_type; x0[ilocal][0] = normal(0); @@ -319,8 +306,7 @@ void FixSMDWallSurface::read_triangles(int pass) { bigint nblocal = atom->nlocal; MPI_Allreduce(&nblocal, &atom->natoms, 1, MPI_LMP_BIGINT, MPI_SUM, world); - if (atom->natoms < 0 || atom->natoms >= MAXBIGINT) - error->all(FLERR, "Too many total atoms"); + if (atom->natoms < 0 || atom->natoms >= MAXBIGINT) error->all(FLERR, "Too many total atoms"); // add IDs for newly created atoms // check that atom IDs are valid @@ -347,3 +333,47 @@ void FixSMDWallSurface::read_triangles(int pass) { delete[] vert; fclose(fp); } + +/* ---------------------------------------------------------------------- + create list of visualization object for dump image + ------------------------------------------------------------------------- */ +int FixSMDWallSurface::image(int *&objs, double **&parms) +{ + const int nlocal = atom->nlocal; + const int * const type = atom->type; + const double * const * const verts = atom->smd_data_9; + + int numobjs = 0; + for (int i = 0; i < nlocal; ++i) + if (type[i] == wall_particle_type) ++numobjs; + + if (numobjs == 0) return 0; + + // reallocate storage + memory->destroy(imgobjs); + memory->destroy(imgparms); + memory->create(imgobjs, numobjs, "wall_surface:imgobjs"); + memory->create(imgparms, numobjs, 10, "wall_surface:imgparms"); + + // copy local tri object info + numobjs = 0; + for (int i = 0; i < nlocal; ++i) { + if (type[i] == wall_particle_type) { + imgobjs[numobjs] = DumpImage::TRI; + imgparms[numobjs][0] = wall_particle_type; + imgparms[numobjs][1] = verts[i][0]; + imgparms[numobjs][2] = verts[i][1]; + imgparms[numobjs][3] = verts[i][2]; + imgparms[numobjs][4] = verts[i][3]; + imgparms[numobjs][5] = verts[i][4]; + imgparms[numobjs][6] = verts[i][5]; + imgparms[numobjs][7] = verts[i][6]; + imgparms[numobjs][8] = verts[i][7]; + imgparms[numobjs][9] = verts[i][8]; + ++numobjs; + } + } + objs = imgobjs; + parms = imgparms; + return numobjs; +} diff --git a/src/MACHDYN/fix_smd_wall_surface.h b/src/MACHDYN/fix_smd_wall_surface.h index fc7478c52bd..609537961f6 100644 --- a/src/MACHDYN/fix_smd_wall_surface.h +++ b/src/MACHDYN/fix_smd_wall_surface.h @@ -37,12 +37,19 @@ class FixSMDWallSurface : public Fix { void read_triangles(int pass); + int image(int *&, double **&) override; + private: int first; // flag for first time initialization double sublo[3], subhi[3]; // epsilon-extended proc sub-box for adding atoms; char *filename; int wall_particle_type; int wall_molecule_id; + + // arrays for dump image rendering + + int *imgobjs; + double **imgparms; }; } // namespace LAMMPS_NS diff --git a/src/REAXFF/fix_reaxff_bonds.cpp b/src/REAXFF/fix_reaxff_bonds.cpp index 74c1b67bbae..f2b1412f48a 100644 --- a/src/REAXFF/fix_reaxff_bonds.cpp +++ b/src/REAXFF/fix_reaxff_bonds.cpp @@ -21,6 +21,7 @@ #include "atom.h" #include "comm.h" #include "error.h" +#include "dump_image.h" #include "force.h" #include "memory.h" #include "neigh_list.h" @@ -40,7 +41,7 @@ using namespace ReaxFF; FixReaxFFBonds::FixReaxFFBonds(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), neighid(nullptr), abo(nullptr), fp(nullptr), lists(nullptr), - reaxff(nullptr), list(nullptr) + reaxff(nullptr), list(nullptr), imgobjs(nullptr), imgparms(nullptr) { if (narg != 5) error->all(FLERR, Error::NOPOINTER, "Fix reaxff/bonds expected 5 arguments but got {}", narg); @@ -50,9 +51,11 @@ FixReaxFFBonds::FixReaxFFBonds(LAMMPS *lmp, int narg, char **arg) : multifile = 0; padflag = 0; first_flag = true; + numobjs = 0; nevery = utils::inumeric(FLERR,arg[3],false,lmp); if (nevery <= 0) error->all(FLERR, 3, "Illegal fix reaxff/bonds nevery value {}", nevery); + global_freq = nevery; filename = arg[4]; if (filename.rfind('*') != std::string::npos) multifile = 1; @@ -71,6 +74,10 @@ FixReaxFFBonds::~FixReaxFFBonds() destroy(); if (fp) fclose(fp); + + // clean up dump image data + memory->destroy(imgobjs); + memory->destroy(imgparms); } /* ---------------------------------------------------------------------- */ @@ -140,6 +147,9 @@ void FixReaxFFBonds::Output_ReaxFF_Bonds() } numbonds = FindBond(); + // no file output with NULL file name. + if (filename == "NULL") return; + // allocate a temporary buffer for the snapshot info MPI_Allreduce(&numbonds,&numbonds_max,1,MPI_INT,MPI_MAX,world); MPI_Allreduce(&nlocal,&nlocal_max,1,MPI_INT,MPI_MAX,world); @@ -155,7 +165,6 @@ void FixReaxFFBonds::Output_ReaxFF_Bonds() RecvBuffer(buf, nbuf, nbuf_local, nlocal_tot, numbonds_max); memory->destroy(buf); - } /* ---------------------------------------------------------------------- */ @@ -174,6 +183,7 @@ int FixReaxFFBonds::FindBond() tagint *tag = atom->tag; int numbonds = 0; + numobjs = 0; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; @@ -192,6 +202,7 @@ int FixReaxFFBonds::FindBond() } } numneigh[i] = nj; + numobjs += nj; if (nj > numbonds) numbonds = nj; } return numbonds; @@ -359,3 +370,49 @@ double FixReaxFFBonds::memory_usage() return bytes; } + +/* ---------------------------------------------------------------------- */ + +int FixReaxFFBonds::image(int *&objs, double **&parms) +{ + if (atom->map_style == Atom::MAP_NONE) + error->all(FLERR, Error::NOLASTLINE, + "Using fix reaxff/bonds with dump image requires an atom map"); + + if (!numobjs) return 0; + memory->destroy(imgobjs); + memory->destroy(imgparms); + memory->create(imgobjs, numobjs, "reaxff/bonds:imgobjs"); + memory->create(imgparms, numobjs, 8, "reaxff/bonds:imgparms"); + + const int nlocal = atom->nlocal; + const int *type = atom->type; + const double * const * const x = atom->x; + + int inum = reaxff->list->inum; + int *ilist = reaxff->list->ilist; + + int n = 0; + for (int ii = 0; ii < inum; ++ii) { + int i = ilist[ii]; + for (int jj = 0; jj < numneigh[i]; ++jj) { + int j = atom->map(neighid[i][jj]); + j = domain->closest_image(i,j); + if (j < 0) continue; + imgobjs[n] = DumpImage::BOND; + imgparms[n][0] = type[i]; + imgparms[n][1] = type[j]; + imgparms[n][2] = x[i][0]; + imgparms[n][3] = x[i][1]; + imgparms[n][4] = x[i][2]; + imgparms[n][5] = x[j][0]; + imgparms[n][6] = x[j][1]; + imgparms[n][7] = x[j][2]; + ++n; + } + } + + objs = imgobjs; + parms = imgparms; + return n; +} diff --git a/src/REAXFF/fix_reaxff_bonds.h b/src/REAXFF/fix_reaxff_bonds.h index 1cc912be391..29cd309b530 100644 --- a/src/REAXFF/fix_reaxff_bonds.h +++ b/src/REAXFF/fix_reaxff_bonds.h @@ -33,6 +33,8 @@ class FixReaxFFBonds : public Fix { void setup(int) override; void end_of_step() override; + int image(int *&, double **&) override; + protected: int nmax, compressed, multifile, padflag; int *numneigh; @@ -54,8 +56,13 @@ class FixReaxFFBonds : public Fix { struct _reax_list *lists; class PairReaxFF *reaxff; class NeighList *list; + + // arrays for dump image rendering + + int numobjs; + int *imgobjs; + double **imgparms; }; } // namespace LAMMPS_NS - #endif #endif diff --git a/src/dump_image.cpp b/src/dump_image.cpp index ad07a2e5d04..86130ca87a7 100644 --- a/src/dump_image.cpp +++ b/src/dump_image.cpp @@ -58,6 +58,8 @@ #include #include +// clang-format on + // helper functions for generating triangle meshes namespace { @@ -204,8 +206,8 @@ void ellipsoid2wireframe(LAMMPS_NS::Image *img, int level, const double *color, } } -void ellipsoid2filled(LAMMPS_NS::Image *img, int level, const double *color, - const double *center, const double *radius, LAMMPS_NS::Region *reg) +void ellipsoid2filled(LAMMPS_NS::Image *img, int level, const double *color, const double *center, + const double *radius, LAMMPS_NS::Region *reg, double opacity) { vec3 offset = {center[0], center[1], center[2]}; @@ -228,7 +230,7 @@ void ellipsoid2filled(LAMMPS_NS::Image *img, int level, const double *color, reg->forward_transform(tri[0][0], tri[0][1], tri[0][2]); reg->forward_transform(tri[1][0], tri[1][1], tri[1][2]); reg->forward_transform(tri[2][0], tri[2][1], tri[2][2]); - img->draw_triangle(tri[0].data(), tri[1].data(), tri[2].data(), color); + img->draw_triangle(tri[0].data(), tri[1].data(), tri[2].data(), color, opacity); } } @@ -242,7 +244,7 @@ void ellipsoid2filled(LAMMPS_NS::Image *img, int level, const double *color, reg->forward_transform(tri[0][0], tri[0][1], tri[0][2]); reg->forward_transform(tri[1][0], tri[1][1], tri[1][2]); reg->forward_transform(tri[2][0], tri[2][1], tri[2][2]); - img->draw_triangle(tri[0].data(), tri[1].data(), tri[2].data(), color); + img->draw_triangle(tri[0].data(), tri[1].data(), tri[2].data(), color, opacity); } } } @@ -256,7 +258,7 @@ void ellipsoid2filled(LAMMPS_NS::Image *img, int level, const double *color, reg->forward_transform(tri[0][0], tri[0][1], tri[0][2]); reg->forward_transform(tri[1][0], tri[1][1], tri[1][2]); reg->forward_transform(tri[2][0], tri[2][1], tri[2][2]); - img->draw_triangle(tri[0].data(), tri[1].data(), tri[2].data(), color); + img->draw_triangle(tri[0].data(), tri[1].data(), tri[2].data(), color, opacity); } } } @@ -269,7 +271,7 @@ void ellipsoid2filled(LAMMPS_NS::Image *img, int level, const double *color, reg->forward_transform(tri[0][0], tri[0][1], tri[0][2]); reg->forward_transform(tri[1][0], tri[1][1], tri[1][2]); reg->forward_transform(tri[2][0], tri[2][1], tri[2][2]); - img->draw_triangle(tri[0].data(), tri[1].data(), tri[2].data(), color); + img->draw_triangle(tri[0].data(), tri[1].data(), tri[2].data(), color, opacity); } } } @@ -280,11 +282,10 @@ using MathConst::DEG2RAD; static constexpr double BIG = 1.0e20; -enum { NUMERIC, ATOM, TYPE, ELEMENT, ATTRIBUTE }; -enum { SPHERE, LINE, TRI }; // also in some Body and Fix child classes +enum { NUMERIC, ATOM, TYPE, ELEMENT, ATTRIBUTE, CONSTANT }; enum { STATIC, DYNAMIC }; enum { NO = 0, YES = 1, AUTO = 2 }; -enum { FILLED, FRAME, POINTS }; +enum { FILLED, FRAME, POINTS, TRANSPARENT }; /* ---------------------------------------------------------------------- */ @@ -292,10 +293,10 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : DumpCustom(lmp, narg, arg), thetastr(nullptr), phistr(nullptr), cxstr(nullptr), cystr(nullptr), czstr(nullptr), upxstr(nullptr), upystr(nullptr), upzstr(nullptr), zoomstr(nullptr), diamtype(nullptr), diamelement(nullptr), bdiamtype(nullptr), colortype(nullptr), - colorelement(nullptr), bcolortype(nullptr), grid2d(nullptr), grid3d(nullptr), - id_grid_compute(nullptr), id_grid_fix(nullptr), grid_compute(nullptr), grid_fix(nullptr), - gbuf(nullptr), avec_line(nullptr), avec_tri(nullptr), avec_body(nullptr), fixptr(nullptr), - image(nullptr), chooseghost(nullptr), bufcopy(nullptr) + colorelement(nullptr), bcolortype(nullptr), aopacity(nullptr), bopacity(nullptr), + grid2d(nullptr), grid3d(nullptr), id_grid_compute(nullptr), id_grid_fix(nullptr), + grid_compute(nullptr), grid_fix(nullptr), gbuf(nullptr), avec_line(nullptr), avec_tri(nullptr), + avec_body(nullptr), image(nullptr), chooseghost(nullptr), bufcopy(nullptr) { if (binary || multiproc) error->all(FLERR, 4, "Invalid dump image filename {}", filename); @@ -308,6 +309,7 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : has_id = true; + // clang-format off // set filetype based on filename suffix if (utils::strmatch(filename, "\\.jpg$") || utils::strmatch(filename, "\\.JPG$") @@ -355,16 +357,16 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : atomflag = YES; gridflag = NO; - lineflag = triflag = bodyflag = fixflag = NO; + lineflag = triflag = bodyflag = NO; - if (atom->nbondtypes == 0) bondflag = NO; - else { + bcolor = ATOM; + bdiam = NUMERIC; + bdiamvalue = 0.5; + if (atom->nbondtypes == 0) { + bondflag = NO; + } else { bondflag = YES; - bcolor = ATOM; - bdiam = NUMERIC; - bdiamvalue = 0.5; } - char *fixID = nullptr; cflag = STATIC; cx = cy = cz = 0.5; @@ -373,6 +375,9 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : boxdiam = 0.02; axesflag = NO; subboxflag = NO; + boxopacity = 1.0; + axesopacity = 1.0; + subboxopacity = 1.0; // parse optional args @@ -477,12 +482,17 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg],"fix") == 0) { if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"dump image fix", error); - fixflag = YES; - fixID = arg[iarg+1]; + std::string id_fix = arg[iarg+1]; + auto *fixptr = modify->get_fix_by_id(id_fix); + if (!fixptr) error->all(FLERR, iarg+1, "Dump image fix ID {} does not exist", id_fix); + int fixcolor = TYPE; if (strcmp(arg[iarg+2],"type") == 0) fixcolor = TYPE; - else error->all(FLERR, iarg+2, "Dump image fix only supports color by type"); - fixflag1 = utils::numeric(FLERR,arg[iarg+3],false,lmp); - fixflag2 = utils::numeric(FLERR,arg[iarg+4],false,lmp); + else if (strcmp(arg[iarg+2],"element") == 0) fixcolor = ELEMENT; + else if (strcmp(arg[iarg+2],"const") == 0) fixcolor = CONSTANT; + else error->all(FLERR, iarg+2, "Unsupported color style for dump image fix {}", arg[iarg+2]); + double fixflag1 = utils::numeric(FLERR,arg[iarg+3],false,lmp); + double fixflag2 = utils::numeric(FLERR,arg[iarg+4],false,lmp); + fixes.emplace_back(id_fix, fixptr, fixcolor, fixflag1, fixflag2, image->color2rgb("red")); iarg += 5; } else if (strcmp(arg[iarg],"region") == 0) { @@ -497,17 +507,19 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[iarg+3],"filled") == 0) drawstyle = FILLED; else if (strcmp(arg[iarg+3],"frame") == 0) drawstyle = FRAME; else if (strcmp(arg[iarg+3],"points") == 0) drawstyle = POINTS; + else if (strcmp(arg[iarg+3],"transparent") == 0) drawstyle = TRANSPARENT; else error->all(FLERR, iarg+3, "Unknown region draw style {}", arg[iarg+3]); double framediam = 0.5; + double opacity = 1.0; int npoints = 0; if (drawstyle == FRAME) { - if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"dump image region", error); + if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"dump image region frame", error); framediam = utils::numeric(FLERR, arg[iarg+4], false, lmp); if (framediam <= 0.0) error->all(FLERR, iarg+4, "Dump image region frame diameter must be > 0.0"); ++iarg; } else if (drawstyle == POINTS) { - if (iarg+6 > narg) utils::missing_cmd_args(FLERR,"dump image region", error); + if (iarg+6 > narg) utils::missing_cmd_args(FLERR,"dump image region points", error); npoints = utils::inumeric(FLERR, arg[iarg+4], false, lmp); if (npoints < 1) error->all(FLERR, iarg+4, "Dump image region number of points must be > 0"); @@ -515,9 +527,15 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : if (framediam <= 0.0) error->all(FLERR, iarg+5, "Dump image region point diameter must be > 0.0"); iarg += 2; + } else if (drawstyle == TRANSPARENT) { + if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"dump image region transparent", error); + opacity = utils::numeric(FLERR, arg[iarg+5], false, lmp); + if ((opacity < 0.0) || (opacity > 1.0)) + error->all(FLERR, iarg+5, "Dump image region opacity must be in the range 0.0 to 1.0"); + iarg += 2; } iarg += 4; - regions.emplace_back(regptr->id, regptr, regcolor, drawstyle, framediam, npoints); + regions.emplace_back(regptr->id, regptr, regcolor, drawstyle, framediam, opacity, npoints); } else if (strcmp(arg[iarg],"size") == 0) { if (iarg+3 > narg) utils::missing_cmd_args(FLERR,"dump image size", error); @@ -666,7 +684,7 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : } else error->all(FLERR,"Unknown dump image keyword {}", arg[iarg]); } - // error checks and setup for lineflag, triflag, bodyflag, fixflag + // error checks and setup for lineflag, triflag, bodyflag if (lineflag) { avec_line = dynamic_cast(atom->style_match("line")); @@ -687,13 +705,6 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : extraflag = 0; if (lineflag || triflag || bodyflag) extraflag = 1; - if (fixflag) { - fixptr = modify->get_fix_by_id(fixID); - if (!fixptr) - error->all(FLERR, Error::NOLASTLINE, "Fix ID {} for dump image does not exist", fixID); - - } - // allocate image buffer now that image size is known image->buffers(); @@ -711,9 +722,11 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : diamelement = new double[ntypes+1]; colortype = new double*[ntypes+1]; colorelement = new double*[ntypes+1]; + aopacity = new double[ntypes+1]; for (int i = 1; i <= ntypes; i++) { diamtype[i] = 1.0; + aopacity[i] = 1.0; if (i % 6 == 1) colortype[i] = image->color2rgb("red"); else if (i % 6 == 2) colortype[i] = image->color2rgb("green"); else if (i % 6 == 3) colortype[i] = image->color2rgb("blue"); @@ -725,8 +738,10 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : if (bondflag == YES) { bdiamtype = new double[atom->nbondtypes+1]; bcolortype = new double*[atom->nbondtypes+1]; + bopacity = new double[atom->nbondtypes+1]; for (int i = 1; i <= atom->nbondtypes; i++) { bdiamtype[i] = 0.5; + bopacity[i] = 1.0; if (i % 6 == 1) bcolortype[i] = image->color2rgb("red"); else if (i % 6 == 2) bcolortype[i] = image->color2rgb("green"); else if (i % 6 == 3) bcolortype[i] = image->color2rgb("blue"); @@ -734,9 +749,6 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : else if (i % 6 == 5) bcolortype[i] = image->color2rgb("aqua"); else if (i % 6 == 0) bcolortype[i] = image->color2rgb("cyan"); } - } else { - bdiamtype = nullptr; - bcolortype = nullptr; } // viewflag = DYNAMIC if any view parameter is dynamic @@ -766,8 +778,10 @@ DumpImage::~DumpImage() delete[] diamelement; delete[] colortype; delete[] colorelement; + delete[] aopacity; delete[] bdiamtype; delete[] bcolortype; + delete[] bopacity; memory->destroy(chooseghost); memory->destroy(bufcopy); memory->destroy(gbuf); @@ -915,6 +929,22 @@ void DumpImage::init_style() ((domain->dimension == 3) && domain->zperiodic && (bondcutoff > domain->zprd))) error->all(FLERR, "Dump image autobond cutoff is larger than periodic domain"); } + + // check if fixes with visualization info still exist + for (auto &ifix : fixes) { + auto *fixptr = modify->get_fix_by_id(ifix.id); + if (!fixptr) + error->all(FLERR, Error::NOLASTLINE, "Fix ID {} for dump image does not exist", ifix.id); + ifix.ptr = fixptr; + + // check if fix data for dump image is available at the required steps. + + int nfreq = fixptr->global_freq; + if ((nfreq == 0) || (nevery % nfreq)) + error->all(FLERR, Error::NOLASTLINE, + "Dump {} and fix {} are not executed at compatible timesteps {}", + style, fixptr->style, utils::errorurl(7)); + } } /* ---------------------------------------------------------------------- */ @@ -1200,9 +1230,9 @@ void DumpImage::create_image() { int i,j,k,m,n,itype,atom1,atom2,imol,iatom,btype,ibonus,drawflag; tagint tagprev; - double diameter,delx,dely,delz; - int *bodyvec,*fixvec; - double **bodyarray,**fixarray; + double diameter,delx,dely,delz,opacity; + int *bodyvec; + double **bodyarray; double *color,*color1,*color2; double *p1,*p2,*p3; double pt1[3],pt2[3],pt3[3]; @@ -1251,7 +1281,7 @@ void DumpImage::create_image() if (bodyflag && body[j] >= 0) drawflag = 0; } - if (drawflag) image->draw_sphere(x[j],color,diameter); + if (drawflag) image->draw_sphere(x[j],color,diameter,aopacity[atom->type[j]]); m += size_one; } @@ -1347,7 +1377,7 @@ void DumpImage::create_image() pt2[1] = x[j][1] - dy; pt2[2] = 0.0; - image->draw_cylinder(pt1,pt2,color,ldiamvalue,3); + image->draw_cylinder(pt1,pt2,color,ldiamvalue,3,aopacity[atom->type[j]]); } } @@ -1371,6 +1401,7 @@ void DumpImage::create_image() if (tcolor == TYPE) { color = colortype[type[j]]; } + double opacity = aopacity[atom->type[j]]; MathExtra::quat_to_mat(avec_tri->bonus[tri[i]].quat,mat); MathExtra::matvec(mat,avec_tri->bonus[tri[i]].c1,pt1); @@ -1380,11 +1411,11 @@ void DumpImage::create_image() MathExtra::add3(pt2,x[i],pt2); MathExtra::add3(pt3,x[i],pt3); - if (tridraw) image->draw_triangle(pt1,pt2,pt3,color); + if (tridraw) image->draw_triangle(pt1,pt2,pt3,color,opacity); if (edgedraw) { - image->draw_cylinder(pt1,pt2,color,tdiamvalue,3); - image->draw_cylinder(pt2,pt3,color,tdiamvalue,3); - image->draw_cylinder(pt3,pt1,color,tdiamvalue,3); + image->draw_cylinder(pt1,pt2,color,tdiamvalue,3,opacity); + image->draw_cylinder(pt2,pt3,color,tdiamvalue,3,opacity); + image->draw_cylinder(pt3,pt1,color,tdiamvalue,3,opacity); } } } @@ -1404,15 +1435,15 @@ void DumpImage::create_image() itype = static_cast(buf[m]); color = colortype[itype]; } + double opacity = aopacity[atom->type[j]]; ibonus = body[j]; n = bptr->image(ibonus,bodyflag1,bodyflag2,bodyvec,bodyarray); for (k = 0; k < n; k++) { if (bodyvec[k] == SPHERE) - image->draw_sphere(bodyarray[k],color,bodyarray[k][3]); + image->draw_sphere(bodyarray[k],color,bodyarray[k][3],opacity); else if (bodyvec[k] == LINE) - image->draw_cylinder(&bodyarray[k][0],&bodyarray[k][3], - color,bodyarray[k][6],3); + image->draw_cylinder(&bodyarray[k][0],&bodyarray[k][3],color,bodyarray[k][6],3,opacity); } m += size_one; @@ -1518,9 +1549,9 @@ void DumpImage::create_image() if (adiam == NUMERIC) { diameter = adiamvalue; } else if (adiam == TYPE) { - diameter = MIN(diamtype[type[atom1]],diamtype[type[atom1]]); + diameter = MIN(diamtype[type[atom1]],diamtype[type[atom2]]); } else if (adiam == ELEMENT) { - diameter = MIN(diamelement[type[atom1]],diamelement[type[atom1]]); + diameter = MIN(diamelement[type[atom1]],diamelement[type[atom2]]); } else if (adiam == ATTRIBUTE) { diameter = MIN(bufcopy[atom1][1],bufcopy[atom2][1]); } @@ -1544,16 +1575,16 @@ void DumpImage::create_image() xmid[1] = x[atom1][1] + 0.5*dely; xmid[2] = x[atom1][2] + 0.5*delz; if (bcolor == ATOM) - image->draw_cylinder(x[atom1],xmid,color1,diameter,3); - else image->draw_cylinder(x[atom1],xmid,color,diameter,3); + image->draw_cylinder(x[atom1],xmid,color1,diameter,3,aopacity[type[atom1]]); + else image->draw_cylinder(x[atom1],xmid,color,diameter,3,bopacity[btype]); xmid[0] = x[atom2][0] - 0.5*delx; xmid[1] = x[atom2][1] - 0.5*dely; xmid[2] = x[atom2][2] - 0.5*delz; if (bcolor == ATOM) - image->draw_cylinder(xmid,x[atom2],color2,diameter,3); - else image->draw_cylinder(xmid,x[atom2],color,diameter,3); + image->draw_cylinder(xmid,x[atom2],color2,diameter,3,aopacity[type[atom1]]); + else image->draw_cylinder(xmid,x[atom2],color,diameter,3,bopacity[btype]); - } else image->draw_cylinder(x[atom1],x[atom2],color,diameter,3); + } else image->draw_cylinder(x[atom1],x[atom2],color,diameter,3,bopacity[btype]); } } } @@ -1666,55 +1697,118 @@ void DumpImage::create_image() xmid[0] = x[atom1][0] + 0.5*dx; xmid[1] = x[atom1][1] + 0.5*dy; xmid[2] = x[atom1][2] + 0.5*dz; - image->draw_cylinder(x[atom1],xmid,color1,diameter,3); + image->draw_cylinder(x[atom1],xmid,color1,diameter,3,aopacity[type[atom1]]); xmid[0] = x[atom2][0] - 0.5*dx; xmid[1] = x[atom2][1] - 0.5*dy; xmid[2] = x[atom2][2] - 0.5*dz; - image->draw_cylinder(xmid,x[atom2],color2,diameter,3); + image->draw_cylinder(xmid,x[atom2],color2,diameter,3,aopacity[type[atom2]]); } } } } } - // render objects provided by a fix + // render objects provided by fixes - if (fixflag) { - int tridraw=0,edgedraw=0; - if (domain->dimension == 3) { - tridraw = 1; - edgedraw = 1; - if ((int) fixflag1 == 2) tridraw = 0; - if ((int) fixflag1 == 1) edgedraw = 0; - } + for (const auto &ifix : fixes) { + int *fixvec = nullptr; + double **fixarray = nullptr; + const int ntypes = atom->ntypes; + n = ifix.ptr->image(fixvec,fixarray); + for (i = 0; i < n; i++) { + if (!fixvec || !fixarray) continue; - n = fixptr->image(fixvec,fixarray); + // set color + if (ifix.colorstyle == TYPE) { + itype = static_cast(fixarray[i][0] - 1.0) % ntypes + 1; + color = colortype[itype]; + } else if (ifix.colorstyle == ELEMENT) { + itype = static_cast(fixarray[i][0] - 1.0) % ntypes + 1; + color = colorelement[itype]; + } else if (ifix.colorstyle == CONSTANT) { + color = ifix.rgb; + } else { + color = image->color2rgb("red"); + } - for (i = 0; i < n; i++) { if (fixvec[i] == SPHERE) { - // no fix draws spheres yet + image->draw_sphere(&fixarray[i][1],color,fixarray[i][4]+ifix.flag2); } else if (fixvec[i] == LINE) { - if (fixcolor == TYPE) { - itype = static_cast(fixarray[i][0]); - color = colortype[itype]; + // @sjplimp for consistency this should be: + // image->draw_cylinder(&fixarray[i][1],&fixarray[i][4],color,ifix.flag2,ifix.flag1); + image->draw_cylinder(&fixarray[i][1],&fixarray[i][4],color,ifix.flag1,3); + } else if (fixvec[i] == TRI) { // don't render surface meshes in 2d + if (domain->dimension == 3) { + p1 = &fixarray[i][1]; + p2 = &fixarray[i][4]; + p3 = &fixarray[i][7]; + if (static_cast(ifix.flag1) % 2) { + image->draw_triangle(p1,p2,p3,color,ifix.flag2); + } else { + image->draw_cylinder(p1,p2,color,ifix.flag2,3); + image->draw_cylinder(p2,p3,color,ifix.flag2,3); + image->draw_cylinder(p3,p1,color,ifix.flag2,3); + } } - image->draw_cylinder(&fixarray[i][1],&fixarray[i][4], - color,fixflag1,3); - } else if (fixvec[i] == TRI) { - if (fixcolor == TYPE) { - itype = static_cast(fixarray[i][0]); - color = colortype[itype]; + } else if (fixvec[i] == CYLINDER) { + image->draw_cylinder(&fixarray[i][1],&fixarray[i][4],color, + fixarray[i][7]+ifix.flag2,(int)ifix.flag1); + } else if (fixvec[i] == TRIANGLE) { + image->draw_triangle(&fixarray[i][1],&fixarray[i][4],&fixarray[i][7],color,ifix.flag2); + } else if (fixvec[i] == BOND) { + int type1 = static_cast(fixarray[i][0] - 1.0) % ntypes + 1; + int type2 = static_cast(fixarray[i][1] - 1.0) % ntypes + 1; + double *color1; + double *color2; + if (ifix.colorstyle == TYPE) { + color1 = colortype[type1]; + color2 = colortype[type2]; + } else if (ifix.colorstyle == ELEMENT) { + color1 = colorelement[type1]; + color2 = colorelement[type2]; + } else if (ifix.colorstyle == CONSTANT) { + color1 = ifix.rgb; + color2 = ifix.rgb; + } else { + color1 = image->color2rgb("white"); + color2 = image->color2rgb("white"); } - p1 = &fixarray[i][1]; - p2 = &fixarray[i][4]; - p3 = &fixarray[i][7]; - if (tridraw) image->draw_triangle(p1,p2,p3,color); - if (edgedraw) { - image->draw_cylinder(p1,p2,color,fixflag2,3); - image->draw_cylinder(p2,p3,color,fixflag2,3); - image->draw_cylinder(p3,p1,color,fixflag2,3); + + double diameter = 0.5; + if (bdiam == ATOM) { + if (adiam == NUMERIC) { + diameter = adiamvalue; + } else if (adiam == TYPE) { + diameter = MIN(diamtype[type1],diamtype[type2]); + } else if (adiam == ELEMENT) { + diameter = MIN(diamelement[type1],diamelement[type2]); + } else if (adiam == ATTRIBUTE) { + diameter = MIN(bufcopy[atom1][1],bufcopy[atom2][1]); + } + } else { + diameter = bdiamvalue; } + // bond diameter adjustment from dump image command line + diameter += ifix.flag2; + + // draw bond cylinder in 2 pieces + + int capflag = ifix.flag1 ? 3 : 0; + delx = fixarray[i][5] - fixarray[i][2]; + dely = fixarray[i][6] - fixarray[i][3]; + delz = fixarray[i][7] - fixarray[i][4]; + + domain->minimum_image(FLERR,delx,dely,delz); + double xmid[3]; + xmid[0] = fixarray[i][2] + 0.5*delx; + xmid[1] = fixarray[i][3] + 0.5*dely; + xmid[2] = fixarray[i][4] + 0.5*delz; + image->draw_cylinder(&fixarray[i][2],xmid,color1,diameter,capflag); + xmid[0] = fixarray[i][5] - 0.5*delx; + xmid[1] = fixarray[i][6] - 0.5*dely; + xmid[2] = fixarray[i][7] - 0.5*delz; + image->draw_cylinder(xmid,&fixarray[i][5],color2,diameter,capflag); } } } @@ -1797,30 +1891,31 @@ void DumpImage::create_image() image->draw_cylinder(block[5],block[6],reg.color,reg.diameter,3); image->draw_cylinder(block[4],block[7],reg.color,reg.diameter,3); image->draw_cylinder(block[6],block[7],reg.color,reg.diameter,3); - } else if (reg.style == FILLED) { + } else if ((reg.style == FILLED) || (reg.style == TRANSPARENT)) { + double opacity = (reg.style == TRANSPARENT) ? reg.opacity : 1.0; if (!myreg->open_faces[0]) { - image->draw_triangle(block[0], block[1], block[2], reg.color); - image->draw_triangle(block[2], block[3], block[0], reg.color); + image->draw_triangle(block[0], block[1], block[2], reg.color, opacity); + image->draw_triangle(block[2], block[3], block[0], reg.color, opacity); } if (!myreg->open_faces[1]) { - image->draw_triangle(block[4], block[5], block[6], reg.color); - image->draw_triangle(block[6], block[7], block[4], reg.color); + image->draw_triangle(block[4], block[5], block[6], reg.color, opacity); + image->draw_triangle(block[6], block[7], block[4], reg.color, opacity); } if (!myreg->open_faces[2]) { - image->draw_triangle(block[0], block[4], block[7], reg.color); - image->draw_triangle(block[7], block[3], block[0], reg.color); + image->draw_triangle(block[0], block[4], block[7], reg.color, opacity); + image->draw_triangle(block[7], block[3], block[0], reg.color, opacity); } if (!myreg->open_faces[3]) { - image->draw_triangle(block[1], block[2], block[6], reg.color); - image->draw_triangle(block[6], block[5], block[1], reg.color); + image->draw_triangle(block[1], block[2], block[6], reg.color, opacity); + image->draw_triangle(block[6], block[5], block[1], reg.color, opacity); } if (!myreg->open_faces[4]) { - image->draw_triangle(block[0], block[1], block[5], reg.color); - image->draw_triangle(block[5], block[4], block[0], reg.color); + image->draw_triangle(block[0], block[1], block[5], reg.color, opacity); + image->draw_triangle(block[5], block[4], block[0], reg.color, opacity); } if (!myreg->open_faces[5]) { - image->draw_triangle(block[3], block[2], block[6], reg.color); - image->draw_triangle(block[6], block[7], block[3], reg.color); + image->draw_triangle(block[3], block[2], block[6], reg.color, opacity); + image->draw_triangle(block[6], block[7], block[3], reg.color, opacity); } } @@ -1919,7 +2014,8 @@ void DumpImage::create_image() image->draw_cylinder(p2, p4, reg.color, reg.diameter, 3); } } - } else if (reg.style == FILLED) { + } else if ((reg.style == FILLED) || (reg.style == TRANSPARENT)) { + double opacity = (reg.style == TRANSPARENT) ? reg.opacity : 1.0; for (int i = 0; i < RESOLUTION; ++i) { if (myreg->axis == 'x') { p1[0] = p2[0] = myreg->lo; @@ -1936,11 +2032,11 @@ void DumpImage::create_image() myreg->forward_transform(p2[0], p2[1], p2[2]); myreg->forward_transform(p3[0], p3[1], p3[2]); myreg->forward_transform(p4[0], p4[1], p4[2]); - if (!myreg->open_faces[0]) image->draw_triangle(p1, p2, lo, reg.color); - if (!myreg->open_faces[1]) image->draw_triangle(p3, p4, hi, reg.color); + if (!myreg->open_faces[0]) image->draw_triangle(p1, p2, lo, reg.color, opacity); + if (!myreg->open_faces[1]) image->draw_triangle(p3, p4, hi, reg.color, opacity); if (!myreg->open_faces[2]) { - image->draw_triangle(p1, p2, p3, reg.color); - image->draw_triangle(p2, p4, p3, reg.color); + image->draw_triangle(p1, p2, p3, reg.color, opacity); + image->draw_triangle(p2, p4, p3, reg.color, opacity); } } else if (myreg->axis == 'y') { @@ -1958,11 +2054,11 @@ void DumpImage::create_image() myreg->forward_transform(p2[0], p2[1], p2[2]); myreg->forward_transform(p3[0], p3[1], p3[2]); myreg->forward_transform(p4[0], p4[1], p4[2]); - if (!myreg->open_faces[0]) image->draw_triangle(p1, p2, lo, reg.color); - if (!myreg->open_faces[1]) image->draw_triangle(p3, p4, hi, reg.color); + if (!myreg->open_faces[0]) image->draw_triangle(p1, p2, lo, reg.color, opacity); + if (!myreg->open_faces[1]) image->draw_triangle(p3, p4, hi, reg.color, opacity); if (!myreg->open_faces[2]) { - image->draw_triangle(p1, p2, p3, reg.color); - image->draw_triangle(p2, p4, p3, reg.color); + image->draw_triangle(p1, p2, p3, reg.color, opacity); + image->draw_triangle(p2, p4, p3, reg.color, opacity); } } else { // if (myreg->axis == 'z') p1[2] = p2[2] = myreg->lo; @@ -1979,11 +2075,11 @@ void DumpImage::create_image() myreg->forward_transform(p2[0], p2[1], p2[2]); myreg->forward_transform(p3[0], p3[1], p3[2]); myreg->forward_transform(p4[0], p4[1], p4[2]); - if (!myreg->open_faces[0]) image->draw_triangle(p1, p2, lo, reg.color); - if (!myreg->open_faces[1]) image->draw_triangle(p3, p4, hi, reg.color); + if (!myreg->open_faces[0]) image->draw_triangle(p1, p2, lo, reg.color, opacity); + if (!myreg->open_faces[1]) image->draw_triangle(p3, p4, hi, reg.color, opacity); if (!myreg->open_faces[2]) { - image->draw_triangle(p1, p2, p3, reg.color); - image->draw_triangle(p2, p4, p3, reg.color); + image->draw_triangle(p1, p2, p3, reg.color, opacity); + image->draw_triangle(p2, p4, p3, reg.color, opacity); } } } @@ -2072,7 +2168,8 @@ void DumpImage::create_image() image->draw_cylinder(p2, p4, reg.color, reg.diameter, 3); } } - } else if (reg.style == FILLED) { + } else if ((reg.style == FILLED) || (reg.style == TRANSPARENT)) { + double opacity = (reg.style == TRANSPARENT) ? reg.opacity : 1.0; for (int i = 0; i < RESOLUTION; ++i) { if (myreg->axis == 'x') { p1[0] = p2[0] = myreg->lo; @@ -2085,11 +2182,11 @@ void DumpImage::create_image() myreg->forward_transform(p2[0], p2[1], p2[2]); myreg->forward_transform(p3[0], p3[1], p3[2]); myreg->forward_transform(p4[0], p4[1], p4[2]); - if (!myreg->open_faces[0]) image->draw_triangle(p1, p2, lo, reg.color); - if (!myreg->open_faces[1]) image->draw_triangle(p3, p4, hi, reg.color); + if (!myreg->open_faces[0]) image->draw_triangle(p1, p2, lo, reg.color, opacity); + if (!myreg->open_faces[1]) image->draw_triangle(p3, p4, hi, reg.color, opacity); if (!myreg->open_faces[2]) { - image->draw_triangle(p1, p2, p3, reg.color); - image->draw_triangle(p3, p4, p2, reg.color); + image->draw_triangle(p1, p2, p3, reg.color, opacity); + image->draw_triangle(p3, p4, p2, reg.color, opacity); } } else if (myreg->axis == 'y') { p1[1] = p2[1] = myreg->lo; @@ -2102,11 +2199,11 @@ void DumpImage::create_image() myreg->forward_transform(p2[0], p2[1], p2[2]); myreg->forward_transform(p3[0], p3[1], p3[2]); myreg->forward_transform(p4[0], p4[1], p4[2]); - if (!myreg->open_faces[0]) image->draw_triangle(p1, p2, lo, reg.color); - if (!myreg->open_faces[1]) image->draw_triangle(p3, p4, hi, reg.color); + if (!myreg->open_faces[0]) image->draw_triangle(p1, p2, lo, reg.color, opacity); + if (!myreg->open_faces[1]) image->draw_triangle(p3, p4, hi, reg.color, opacity); if (!myreg->open_faces[2]) { - image->draw_triangle(p1, p2, p3, reg.color); - image->draw_triangle(p3, p4, p2, reg.color); + image->draw_triangle(p1, p2, p3, reg.color, opacity); + image->draw_triangle(p3, p4, p2, reg.color, opacity); } } else { // if (myreg->axis == 'z') p1[2] = p2[2] = myreg->lo; @@ -2119,11 +2216,11 @@ void DumpImage::create_image() myreg->forward_transform(p2[0], p2[1], p2[2]); myreg->forward_transform(p3[0], p3[1], p3[2]); myreg->forward_transform(p4[0], p4[1], p4[2]); - if (!myreg->open_faces[0]) image->draw_triangle(p1, p2, lo, reg.color); - if (!myreg->open_faces[1]) image->draw_triangle(p3, p4, hi, reg.color); + if (!myreg->open_faces[0]) image->draw_triangle(p1, p2, lo, reg.color, opacity); + if (!myreg->open_faces[1]) image->draw_triangle(p3, p4, hi, reg.color, opacity); if (!myreg->open_faces[2]) { - image->draw_triangle(p1, p2, p3, reg.color); - image->draw_triangle(p2, p4, p3, reg.color); + image->draw_triangle(p1, p2, p3, reg.color, opacity); + image->draw_triangle(p2, p4, p3, reg.color, opacity); } } } @@ -2141,8 +2238,9 @@ void DumpImage::create_image() double radius[3] = {myreg->a, myreg->b, myreg->c}; if (reg.style == FRAME) { ellipsoid2wireframe(image, 4, reg.color, reg.diameter, center, radius, reg.ptr); - } else if (reg.style == FILLED) { - ellipsoid2filled(image, 4, reg.color, center, radius, reg.ptr); + } else if ((reg.style == FILLED) || (reg.style == TRANSPARENT)) { + double opacity = (reg.style == TRANSPARENT) ? reg.opacity : 1.0; + ellipsoid2filled(image, 4, reg.color, center, radius, reg.ptr, opacity); } } else if (regstyle == "prism") { @@ -2175,30 +2273,31 @@ void DumpImage::create_image() image->draw_cylinder(block[5],block[6],reg.color,reg.diameter,3); image->draw_cylinder(block[4],block[7],reg.color,reg.diameter,3); image->draw_cylinder(block[6],block[7],reg.color,reg.diameter,3); - } else if (reg.style == FILLED) { + } else if ((reg.style == FILLED) || (reg.style == TRANSPARENT)) { + double opacity = (reg.style == TRANSPARENT) ? reg.opacity : 1.0; if (!myreg->open_faces[0]) { - image->draw_triangle(block[0], block[1], block[2], reg.color); - image->draw_triangle(block[2], block[3], block[0], reg.color); + image->draw_triangle(block[0], block[1], block[2], reg.color, opacity); + image->draw_triangle(block[2], block[3], block[0], reg.color, opacity); } if (!myreg->open_faces[1]) { - image->draw_triangle(block[4], block[5], block[6], reg.color); - image->draw_triangle(block[6], block[7], block[4], reg.color); + image->draw_triangle(block[4], block[5], block[6], reg.color, opacity); + image->draw_triangle(block[6], block[7], block[4], reg.color, opacity); } if (!myreg->open_faces[2]) { - image->draw_triangle(block[0], block[4], block[7], reg.color); - image->draw_triangle(block[7], block[3], block[0], reg.color); + image->draw_triangle(block[0], block[4], block[7], reg.color, opacity); + image->draw_triangle(block[7], block[3], block[0], reg.color, opacity); } if (!myreg->open_faces[3]) { - image->draw_triangle(block[1], block[2], block[6], reg.color); - image->draw_triangle(block[6], block[5], block[1], reg.color); + image->draw_triangle(block[1], block[2], block[6], reg.color, opacity); + image->draw_triangle(block[6], block[5], block[1], reg.color, opacity); } if (!myreg->open_faces[4]) { - image->draw_triangle(block[0], block[1], block[5], reg.color); - image->draw_triangle(block[5], block[4], block[0], reg.color); + image->draw_triangle(block[0], block[1], block[5], reg.color, opacity); + image->draw_triangle(block[5], block[4], block[0], reg.color, opacity); } if (!myreg->open_faces[5]) { - image->draw_triangle(block[3], block[2], block[6], reg.color); - image->draw_triangle(block[6], block[7], block[3], reg.color); + image->draw_triangle(block[3], block[2], block[6], reg.color, opacity); + image->draw_triangle(block[6], block[7], block[3], reg.color, opacity); } } } else if (regstyle == "sphere") { @@ -2213,9 +2312,10 @@ void DumpImage::create_image() if (reg.style == FRAME) { double radius[3] = {myreg->radius,myreg->radius,myreg->radius}; ellipsoid2wireframe(image, 4, reg.color, reg.diameter, center, radius, reg.ptr); - } else if (reg.style == FILLED) { + } else if ((reg.style == FILLED) || (reg.style == TRANSPARENT)) { + double opacity = (reg.style == TRANSPARENT) ? reg.opacity : 1.0; myreg->forward_transform(center[0], center[1], center[2]); - image->draw_sphere(center, reg.color, 2.0 * myreg->radius); + image->draw_sphere(center, reg.color, 2.0 * myreg->radius, opacity); } } else { if (comm->me == 0) @@ -2251,7 +2351,7 @@ void DumpImage::create_image() boxcorners = domain->corners; } - image->draw_box(boxcorners,diameter); + image->draw_box(boxcorners,diameter,subboxopacity); } // render outline of simulation box, orthogonal or triclinic @@ -2278,7 +2378,7 @@ void DumpImage::create_image() boxcorners = domain->corners; } - image->draw_box(boxcorners,diameter); + image->draw_box(boxcorners,diameter,boxopacity); } // render XYZ axes in red/green/blue @@ -2331,7 +2431,7 @@ void DumpImage::create_image() axes[3][1] = axes[0][1] + axeslen*(axes[3][1]-axes[0][1]); axes[3][2] = axes[0][2] + axeslen*(axes[3][2]-axes[0][2]); - image->draw_axes(axes,diameter); + image->draw_axes(axes,diameter,axesopacity); } } @@ -2499,6 +2599,16 @@ int DumpImage::modify_param(int narg, char **arg) return 3; } + if (strcmp(arg[0],"atrans") == 0) { + if (narg < 3) error->all(FLERR,"Illegal dump_modify command"); + int nlo,nhi; + utils::bounds_typelabel(FLERR,arg[1],1,atom->ntypes,nlo,nhi,lmp,Atom::ATOM); + double opacity = utils::numeric(FLERR,arg[2],false,lmp); + if ((opacity < 0.0) || (opacity > 1.0)) error->all(FLERR,"Illegal dump_modify command"); + for (int i = nlo; i <= nhi; i++) aopacity[i] = opacity; + return 3; + } + if ((strcmp(arg[0],"amap") == 0) || (strcmp(arg[0],"gmap") == 0)) { if (narg < 6) error->all(FLERR,"Illegal dump_modify command"); if (strlen(arg[3]) != 2) error->all(FLERR,"Illegal dump_modify command"); @@ -2553,6 +2663,18 @@ int DumpImage::modify_param(int narg, char **arg) return 3; } + if (strcmp(arg[0],"btrans") == 0) { + if (narg < 3) error->all(FLERR,"Illegal dump_modify command"); + if (atom->nbondtypes == 0) + error->all(FLERR,"Dump modify btrans not allowed with no bond types"); + int nlo,nhi; + utils::bounds_typelabel(FLERR,arg[1],1,atom->nbondtypes,nlo,nhi,lmp,Atom::BOND); + double opacity = utils::numeric(FLERR,arg[2],false,lmp); + if ((opacity < 0.0) || (opacity > 1.0)) error->all(FLERR,"Illegal dump_modify command"); + for (int i = nlo; i <= nhi; i++) bopacity[i] = opacity; + return 3; + } + if (strcmp(arg[0],"backcolor") == 0) { if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); double *color = image->color2rgb(arg[1]); @@ -2571,6 +2693,30 @@ int DumpImage::modify_param(int narg, char **arg) return 2; } + if (strcmp(arg[0],"boxtrans") == 0) { + if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); + boxopacity = utils::numeric(FLERR,arg[1],false,lmp); + if ((boxopacity < 0.0) || (boxopacity > 1.0)) + error->all(FLERR,"Invalid boxtrans in dump_modify command"); + return 2; + } + + if (strcmp(arg[0],"subboxtrans") == 0) { + if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); + subboxopacity = utils::numeric(FLERR,arg[1],false,lmp); + if ((subboxopacity < 0.0) || (subboxopacity > 1.0)) + error->all(FLERR,"Invalid subboxtrans in dump_modify command"); + return 2; + } + + if (strcmp(arg[0],"axestrans") == 0) { + if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); + axesopacity = utils::numeric(FLERR,arg[1],false,lmp); + if ((axesopacity < 0.0) || (axesopacity > 1.0)) + error->all(FLERR,"Invalid axestrans in dump_modify command"); + return 2; + } + if (strcmp(arg[0],"color") == 0) { if (narg < 5) error->all(FLERR,"Illegal dump_modify command"); int flag = image->addcolor(arg[1],utils::numeric(FLERR,arg[2],false,lmp), @@ -2580,5 +2726,15 @@ int DumpImage::modify_param(int narg, char **arg) return 5; } + if (strcmp(arg[0],"fcolor") == 0) { + if (narg < 3) error->all(FLERR,"Illegal dump_modify command"); + auto *color = image->color2rgb(arg[2]); + if (!color) error->all(FLERR, "Unknown color for dump_modify fcolor: {}", arg[2]); + for (auto &ifix : fixes) { + if (ifix.id == arg[1]) ifix.rgb = color; + } + return 3; + } + return 0; } diff --git a/src/dump_image.h b/src/dump_image.h index 1cb13c6ff2d..80efc057543 100644 --- a/src/dump_image.h +++ b/src/dump_image.h @@ -23,19 +23,39 @@ DumpStyle(image,DumpImage); #include "dump_custom.h" namespace LAMMPS_NS { + +// forward declarations +class AtomVecBody; +class AtomVecLine; +class AtomVecTri; +class Compute; +class Fix; +class Grid2d; +class Grid3d; +class Image; +class LAMMPS; class Region; + class DumpImage : public DumpCustom { public: int multifile_override; // used by write_dump command - - DumpImage(class LAMMPS *, int, char **); + enum { + SPHERE, // a single sphere with radius provided + LINE, // a cylinder with diameter given through fflag2 + TRI, // a surface mesh as triangles or cylinder mesh based on fflag1, fflag2 sets diameter + CYLINDER, // a cylinder with diameter given by fix, fflag1 choose caps, fflag2 adjusts diameter + TRIANGLE, // a regular triangle, no settings apply + BOND // two connected cylinders with bond diameter, colored by atom types, fflag1 sets cap + }; // used by some Body and Fix child classes + + DumpImage(LAMMPS *, int, char **); ~DumpImage() override; int pack_forward_comm(int, int *, double *, int, int *) override; void unpack_forward_comm(int, int, double *) override; protected: int filetype; - enum { PPM, JPG, PNG }; + enum { PPM, JPG, PNG }; // file type constants int atomflag; // 0/1 for draw atoms int acolor, adiam; // what determines color/diam of atoms @@ -50,9 +70,6 @@ class DumpImage : public DumpCustom { int bodyflag; // 0/1 for draw atoms as bodies int bodycolor; // what determines color of bodies double bodyflag1, bodyflag2; // user-specified params for drawing bodies - int fixflag; // 0/1 to draw what fix provides - int fixcolor; // what determines color of fix objects - double fixflag1, fixflag2; // user-specified params for fix objects int bondflag; // NO/YES/AUTO for drawing bonds int bcolor, bdiam; // what determines color/diam of bonds @@ -72,6 +89,7 @@ class DumpImage : public DumpCustom { int zoomvar; // index to zoom variable int boxflag, axesflag; // 0/1 for draw box and axes double boxdiam, axeslen, axesdiam; // params for drawing box and axes + double boxopacity, axesopacity, subboxopacity; // opacity for box, subbox, axes int subboxflag; double subboxdiam; @@ -79,13 +97,14 @@ class DumpImage : public DumpCustom { double *diamtype, *diamelement, *bdiamtype; // per-type diameters double **colortype, **colorelement, **bcolortype; // per-type colors + double *aopacity, *bopacity; // per-type opacity int gridflag; // 0/1 for draw grid cells - class Grid2d *grid2d; - class Grid3d *grid3d; + Grid2d *grid2d; + Grid3d *grid3d; char *id_grid_compute, *id_grid_fix; - class Compute *grid_compute; - class Fix *grid_fix; + Compute *grid_compute; + Fix *grid_fix; int grid_igrid, grid_idata, grid_index; int nxgrid, nygrid, nzgrid; int nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in; @@ -93,20 +112,36 @@ class DumpImage : public DumpCustom { int ngrid, maxgrid; double gcorners[8][3]; - class AtomVecLine *avec_line; // ptrs to atom style (sub)classes - class AtomVecTri *avec_tri; - class AtomVecBody *avec_body; + AtomVecLine *avec_line; // ptrs to atom style (sub)classes + AtomVecTri *avec_tri; + AtomVecBody *avec_body; + + struct FixInfo { + FixInfo() = delete; + FixInfo(const std::string &_id, Fix *_ptr, int _colorstyle, double _flag1, double _flag2, + double *_rgb) : + id(_id), ptr(_ptr), colorstyle(_colorstyle), flag1(_flag1), flag2(_flag2), rgb(_rgb) + { + } + + Fix *ptr; + std::string id; + int colorstyle; + double flag1; + double flag2; + double *rgb; + }; - class Fix *fixptr; // ptr to Fix that provides image data + std::vector fixes; - class Image *image; // class that renders each image + Image *image; // class that renders each image - class RegionInfo { - public: + struct RegionInfo { RegionInfo() = delete; RegionInfo(const std::string &_id, Region *_ptr, double *_color, int _style, - double _diameter = 0.5, int _npoints = 0) : - ptr(_ptr), id(_id), style(_style), color(_color), diameter(_diameter), npoints(_npoints) + double _diameter = 0.5, double _opacity = 1.0, int _npoints = 0) : + ptr(_ptr), id(_id), style(_style), color(_color), diameter(_diameter), opacity(_opacity), + npoints(_npoints) { } @@ -115,6 +150,7 @@ class DumpImage : public DumpCustom { int style; double *color; double diameter; + double opacity; int npoints; }; diff --git a/src/fix.cpp b/src/fix.cpp index c91990897f9..88d467e8d8e 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -84,7 +84,7 @@ Fix::Fix(LAMMPS *lmp, int /*narg*/, char **arg) : scalar_flag = vector_flag = array_flag = 0; extscalar = extvector = extarray = -1; peratom_flag = local_flag = pergrid_flag = 0; - global_freq = local_freq = peratom_freq = pergrid_freq = -1; + local_freq = peratom_freq = pergrid_freq = -1; size_vector_variable = size_array_rows_variable = 0; comm_forward = comm_reverse = comm_border = 0; diff --git a/src/fix_graphics.cpp b/src/fix_graphics.cpp new file mode 100644 index 00000000000..a58abae62f8 --- /dev/null +++ b/src/fix_graphics.cpp @@ -0,0 +1,680 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "fix_graphics.h" + +#include "comm.h" +#include "domain.h" +#include "dump_image.h" +#include "error.h" +#include "input.h" +#include "lattice.h" +#include "math_extra.h" +#include "memory.h" +#include "modify.h" +#include "update.h" +#include "variable.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum { SPHERE, CYLINDER, ARROW, PROGBAR }; +enum { X = 0, Y, Z }; + +/* ---------------------------------------------------------------------- */ + +FixGraphics::FixGraphics(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), imgobjs(nullptr), imgparms(nullptr) +{ + if (narg < 4) utils::missing_cmd_args(FLERR, "fix graphics", error); + + // parse mandatory arg + + nevery = utils::inumeric(FLERR, arg[3], false, lmp); + if (nevery <= 0) error->all(FLERR, 3, "Illegal fix graphics nevery value {}", nevery); + global_freq = nevery; + dynamic_group_allow = 1; + + numobjs = 0; + varflag = 0; + + int iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg], "sphere") == 0) { + if (iarg + 6 > narg) utils::missing_cmd_args(FLERR, "fix graphics sphere", error); + // clang-format off + SphereItem sphere{SPHERE, 1, {0.0, 0.0, 0.0}, 0.0, nullptr, nullptr, nullptr, nullptr, + -1, -1, -1, -1}; + // clang-format on + sphere.type = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + if (strstr(arg[iarg + 2], "v_") == arg[iarg + 2]) { + varflag = 1; + sphere.xstr = utils::strdup(arg[iarg + 2] + 2); + } else + sphere.pos[X] = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + if (strstr(arg[iarg + 3], "v_") == arg[iarg + 3]) { + varflag = 1; + sphere.ystr = utils::strdup(arg[iarg + 3] + 2); + } else + sphere.pos[Y] = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + if (strstr(arg[iarg + 4], "v_") == arg[iarg + 4]) { + varflag = 1; + sphere.zstr = utils::strdup(arg[iarg + 4] + 2); + } else + sphere.pos[Z] = utils::numeric(FLERR, arg[iarg + 4], false, lmp); + if (strstr(arg[iarg + 5], "v_") == arg[iarg + 5]) { + varflag = 1; + sphere.dstr = utils::strdup(arg[iarg + 5] + 2); + } else + sphere.diameter = 2.0 * utils::numeric(FLERR, arg[iarg + 5], false, lmp); + items.emplace_back(std::move(sphere)); + ++numobjs; + iarg += 6; + } else if (strcmp(arg[iarg], "cylinder") == 0) { + if (iarg + 9 > narg) utils::missing_cmd_args(FLERR, "fix graphics cylinder", error); + // clang-format off + CylinderItem cylinder{CYLINDER, 1, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, 0.0, + nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, + -1, -1, -1, -1, -1, -1, -1}; + // clang-format on + cylinder.type = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + if (strstr(arg[iarg + 2], "v_") == arg[iarg + 2]) { + varflag = 1; + cylinder.x1str = utils::strdup(arg[iarg + 2] + 2); + } else + cylinder.pos1[X] = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + if (strstr(arg[iarg + 3], "v_") == arg[iarg + 3]) { + varflag = 1; + cylinder.y1str = utils::strdup(arg[iarg + 3] + 2); + } else + cylinder.pos1[Y] = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + if (strstr(arg[iarg + 4], "v_") == arg[iarg + 4]) { + varflag = 1; + cylinder.z1str = utils::strdup(arg[iarg + 4] + 2); + } else + cylinder.pos1[Z] = utils::numeric(FLERR, arg[iarg + 4], false, lmp); + if (strstr(arg[iarg + 5], "v_") == arg[iarg + 5]) { + varflag = 1; + cylinder.x2str = utils::strdup(arg[iarg + 5] + 2); + } else + cylinder.pos2[X] = utils::numeric(FLERR, arg[iarg + 5], false, lmp); + if (strstr(arg[iarg + 6], "v_") == arg[iarg + 6]) { + varflag = 1; + cylinder.y2str = utils::strdup(arg[iarg + 6] + 2); + } else + cylinder.pos2[Y] = utils::numeric(FLERR, arg[iarg + 6], false, lmp); + if (strstr(arg[iarg + 7], "v_") == arg[iarg + 7]) { + varflag = 1; + cylinder.z2str = utils::strdup(arg[iarg + 7] + 2); + } else + cylinder.pos2[Z] = utils::numeric(FLERR, arg[iarg + 7], false, lmp); + if (strstr(arg[iarg + 8], "v_") == arg[iarg + 8]) { + varflag = 1; + cylinder.dstr = utils::strdup(arg[iarg + 8] + 2); + } else + cylinder.diameter = 2.0 * utils::numeric(FLERR, arg[iarg + 8], false, lmp); + items.emplace_back(std::move(cylinder)); + ++numobjs; + iarg += 9; + } else if (strcmp(arg[iarg], "arrow") == 0) { + if (iarg + 10 > narg) utils::missing_cmd_args(FLERR, "fix graphics arrow", error); + // clang-format off + ArrowItem arrow{ARROW, 1, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, 0.0, 0.1, + nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, + -1, -1, -1, -1, -1, -1, -1}; + // clang-format on + arrow.type = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + if (strstr(arg[iarg + 2], "v_") == arg[iarg + 2]) { + varflag = 1; + arrow.x1str = utils::strdup(arg[iarg + 2] + 2); + } else + arrow.tip[X] = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + if (strstr(arg[iarg + 3], "v_") == arg[iarg + 3]) { + varflag = 1; + arrow.y1str = utils::strdup(arg[iarg + 3] + 2); + } else + arrow.tip[Y] = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + if (strstr(arg[iarg + 4], "v_") == arg[iarg + 4]) { + varflag = 1; + arrow.z1str = utils::strdup(arg[iarg + 4] + 2); + } else + arrow.tip[Z] = utils::numeric(FLERR, arg[iarg + 4], false, lmp); + if (strstr(arg[iarg + 5], "v_") == arg[iarg + 5]) { + varflag = 1; + arrow.x2str = utils::strdup(arg[iarg + 5] + 2); + } else + arrow.bot[X] = utils::numeric(FLERR, arg[iarg + 5], false, lmp); + if (strstr(arg[iarg + 6], "v_") == arg[iarg + 6]) { + varflag = 1; + arrow.y2str = utils::strdup(arg[iarg + 6] + 2); + } else + arrow.bot[Y] = utils::numeric(FLERR, arg[iarg + 6], false, lmp); + if (strstr(arg[iarg + 7], "v_") == arg[iarg + 7]) { + varflag = 1; + arrow.z2str = utils::strdup(arg[iarg + 7] + 2); + } else + arrow.bot[Z] = utils::numeric(FLERR, arg[iarg + 7], false, lmp); + if (strstr(arg[iarg + 8], "v_") == arg[iarg + 8]) { + varflag = 1; + arrow.dstr = utils::strdup(arg[iarg + 8] + 2); + } else + arrow.diameter = 2.0 * utils::numeric(FLERR, arg[iarg + 8], false, lmp); + arrow.ratio = utils::numeric(FLERR, arg[iarg + 9], false, lmp); + if ((arrow.ratio < 0.1) || (arrow.ratio > 0.5)) + error->all(FLERR, iarg + 9, "Arrow tip ratio must be between 0.1 and 0.5"); + items.emplace_back(std::move(arrow)); + numobjs += 2; + iarg += 10; + } else if (strcmp(arg[iarg], "progbar") == 0) { + if (iarg + 11 > narg) utils::missing_cmd_args(FLERR, "fix graphics progbar", error); + // clang-format off + ProgbarItem progbar{PROGBAR, 1, 2, Y, 0, {0.0, 0.0, 0.0}, 0.0, 0.0, 0.0, nullptr, -1}; + // clang-format on + progbar.type1 = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + progbar.type2 = utils::inumeric(FLERR, arg[iarg + 2], false, lmp); + if (strcmp(arg[iarg + 3], "x") == 0) { + progbar.dim = X; + } else if (strcmp(arg[iarg + 3], "y") == 0) { + progbar.dim = Y; + } else if (strcmp(arg[iarg + 3], "z") == 0) { + progbar.dim = Z; + } else { + error->all(FLERR, iarg + 3, "Unsupported progress bar dimension string {}", arg[iarg + 3]); + } + progbar.pos[X] = utils::numeric(FLERR, arg[iarg + 4], false, lmp); + progbar.pos[Y] = utils::numeric(FLERR, arg[iarg + 5], false, lmp); + progbar.pos[Z] = utils::numeric(FLERR, arg[iarg + 6], false, lmp); + progbar.length = utils::numeric(FLERR, arg[iarg + 7], false, lmp); + if ((progbar.length <= 0.0) || (progbar.length > 2.0 * domain->prd[progbar.dim])) + error->all(FLERR, iarg + 7, "Illegal progress bar length {}", arg[iarg + 7]); + progbar.diameter = 2.0 * utils::numeric(FLERR, arg[iarg + 8], false, lmp); + if (strstr(arg[iarg + 9], "v_") == arg[iarg + 9]) { + varflag = 1; + progbar.pstr = utils::strdup(arg[iarg + 9] + 2); + } else { + progbar.progress = utils::numeric(FLERR, arg[iarg + 9], false, lmp); + } + progbar.tics = utils::inumeric(FLERR, arg[iarg + 10], false, lmp); + if ((progbar.tics < 0) || (progbar.tics > 20)) + error->all(FLERR, iarg + 10, "Unsupported number of progress bar tics {}", arg[iarg + 10]); + items.emplace_back(std::move(progbar)); + numobjs += 2 + progbar.tics; + iarg += 11; + } else { + error->all(FLERR, iarg, "Unknown fix graphics keyword {}", arg[iarg]); + } + } + memory->create(imgobjs, numobjs, "fix_graphics:imgobjs"); + memory->create(imgparms, numobjs, 8, "fix_graphics:imgparms"); +} + +/* ---------------------------------------------------------------------- */ + +FixGraphics::~FixGraphics() +{ + for (auto &gi : items) { + switch (gi.style) { + case SPHERE: + delete[] gi.sphere.xstr; + delete[] gi.sphere.ystr; + delete[] gi.sphere.zstr; + delete[] gi.sphere.dstr; + break; + case CYLINDER: + delete[] gi.cylinder.x1str; + delete[] gi.cylinder.y1str; + delete[] gi.cylinder.z1str; + delete[] gi.cylinder.x2str; + delete[] gi.cylinder.y2str; + delete[] gi.cylinder.z2str; + delete[] gi.cylinder.dstr; + break; + case ARROW: + delete[] gi.arrow.x1str; + delete[] gi.arrow.y1str; + delete[] gi.arrow.z1str; + delete[] gi.arrow.x2str; + delete[] gi.arrow.y2str; + delete[] gi.arrow.z2str; + delete[] gi.arrow.dstr; + break; + case PROGBAR: + delete[] gi.progbar.pstr; + break; + default:; // do nothing + break; + } + } + + memory->destroy(imgobjs); + memory->destroy(imgparms); +} + +/* ---------------------------------------------------------------------- */ + +int FixGraphics::setmask() +{ + return END_OF_STEP; +} + +/* ---------------------------------------------------------------------- */ + +void FixGraphics::init() +{ + int n = 0; + for (auto &gi : items) { + if (gi.style == SPHERE) { + imgobjs[n] = DumpImage::SPHERE; + imgparms[n][0] = gi.sphere.type; + if (gi.sphere.xstr) { + int ivar = input->variable->find(gi.sphere.xstr); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.sphere.xstr); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.sphere.xstr); + gi.sphere.xvar = ivar; + } + if (gi.sphere.ystr) { + int ivar = input->variable->find(gi.sphere.ystr); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.sphere.ystr); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.sphere.ystr); + gi.sphere.yvar = ivar; + } + if (gi.sphere.zstr) { + int ivar = input->variable->find(gi.sphere.zstr); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.sphere.zstr); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.sphere.zstr); + gi.sphere.zvar = ivar; + } + if (gi.sphere.dstr) { + int ivar = input->variable->find(gi.sphere.dstr); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.sphere.dstr); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.sphere.dstr); + gi.sphere.dvar = ivar; + } + ++n; + } else if (gi.style == CYLINDER) { + imgobjs[n] = DumpImage::CYLINDER; + imgparms[n][0] = gi.cylinder.type; + if (gi.cylinder.x1str) { + int ivar = input->variable->find(gi.cylinder.x1str); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.cylinder.x1str); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.cylinder.x1str); + gi.cylinder.x1var = ivar; + } + if (gi.cylinder.y1str) { + int ivar = input->variable->find(gi.cylinder.y1str); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.cylinder.y1str); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.cylinder.y1str); + gi.cylinder.y1var = ivar; + } + if (gi.cylinder.z1str) { + int ivar = input->variable->find(gi.cylinder.z1str); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.cylinder.z1str); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.cylinder.z1str); + gi.cylinder.z1var = ivar; + } + if (gi.cylinder.x2str) { + int ivar = input->variable->find(gi.cylinder.x2str); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.cylinder.x2str); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.cylinder.x2str); + gi.cylinder.x2var = ivar; + } + if (gi.cylinder.y2str) { + int ivar = input->variable->find(gi.cylinder.y2str); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.cylinder.y2str); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.cylinder.y2str); + gi.cylinder.y2var = ivar; + } + if (gi.cylinder.z2str) { + int ivar = input->variable->find(gi.cylinder.z2str); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.cylinder.z2str); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.cylinder.z2str); + gi.cylinder.z2var = ivar; + } + if (gi.cylinder.dstr) { + int ivar = input->variable->find(gi.cylinder.dstr); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.cylinder.dstr); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.cylinder.dstr); + gi.cylinder.dvar = ivar; + } + ++n; + } else if (gi.style == ARROW) { + imgobjs[n] = DumpImage::CYLINDER; + imgparms[n][0] = gi.arrow.type; + if (gi.arrow.x1str) { + int ivar = input->variable->find(gi.arrow.x1str); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.arrow.x1str); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.arrow.x1str); + gi.arrow.x1var = ivar; + } + if (gi.arrow.y1str) { + int ivar = input->variable->find(gi.arrow.y1str); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.arrow.y1str); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.arrow.y1str); + gi.arrow.y1var = ivar; + } + if (gi.arrow.z1str) { + int ivar = input->variable->find(gi.arrow.z1str); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.arrow.z1str); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.arrow.z1str); + gi.arrow.z1var = ivar; + } + if (gi.arrow.x2str) { + int ivar = input->variable->find(gi.arrow.x2str); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.arrow.x2str); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.arrow.x2str); + gi.arrow.x2var = ivar; + } + if (gi.arrow.y2str) { + int ivar = input->variable->find(gi.arrow.y2str); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.arrow.y2str); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.arrow.y2str); + gi.arrow.y2var = ivar; + } + if (gi.arrow.z2str) { + int ivar = input->variable->find(gi.arrow.z2str); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.arrow.z2str); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.arrow.z2str); + gi.arrow.z2var = ivar; + } + if (gi.arrow.dstr) { + int ivar = input->variable->find(gi.arrow.dstr); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.arrow.dstr); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.arrow.dstr); + gi.arrow.dvar = ivar; + } + ++n; + imgobjs[n] = DumpImage::CYLINDER; + imgparms[n][0] = gi.arrow.type; + ++n; + } else if (gi.style == PROGBAR) { + imgobjs[n] = DumpImage::CYLINDER; + imgparms[n][0] = gi.progbar.type1; + imgparms[n][1] = gi.progbar.pos[X]; + imgparms[n][2] = gi.progbar.pos[Y]; + imgparms[n][3] = gi.progbar.pos[Z]; + imgparms[n][4] = gi.progbar.pos[X]; + imgparms[n][5] = gi.progbar.pos[Y]; + imgparms[n][6] = gi.progbar.pos[Z]; + imgparms[n][7] = gi.progbar.diameter; + switch (gi.progbar.dim) { + case X: + imgparms[n][1] -= 0.5 * gi.progbar.length; + imgparms[n][4] += 0.5 * gi.progbar.length; + break; + case Y: + imgparms[n][2] -= 0.5 * gi.progbar.length; + imgparms[n][5] += 0.5 * gi.progbar.length; + break; + case Z: + imgparms[n][3] -= 0.5 * gi.progbar.length; + imgparms[n][6] += 0.5 * gi.progbar.length; + break; + default:; // do nothing + } + ++n; + imgobjs[n] = DumpImage::CYLINDER; + imgparms[n][0] = gi.progbar.type2; + imgparms[n][1] = gi.progbar.pos[X]; + imgparms[n][2] = gi.progbar.pos[Y]; + imgparms[n][3] = gi.progbar.pos[Z]; + imgparms[n][4] = gi.progbar.pos[X]; + imgparms[n][5] = gi.progbar.pos[Y]; + imgparms[n][6] = gi.progbar.pos[Z]; + imgparms[n][7] = 0.75 * gi.progbar.diameter; + switch (gi.progbar.dim) { + case X: + imgparms[n][1] -= 0.5 * gi.progbar.length; + imgparms[n][4] -= 0.5 * gi.progbar.length; + imgparms[n][3] += 0.2 * gi.progbar.diameter; + imgparms[n][6] += 0.2 * gi.progbar.diameter; + break; + case Y: + imgparms[n][2] -= 0.5 * gi.progbar.length; + imgparms[n][5] -= 0.5 * gi.progbar.length; + imgparms[n][1] += 0.15 * gi.progbar.diameter; + imgparms[n][4] += 0.15 * gi.progbar.diameter; + break; + case Z: + imgparms[n][3] -= 0.5 * gi.progbar.length; + imgparms[n][6] -= 0.5 * gi.progbar.length; + imgparms[n][1] += 0.15 * gi.progbar.diameter; + imgparms[n][4] += 0.15 * gi.progbar.diameter; + break; + default: + break; + } + ++n; + double delta = gi.progbar.length / (double) (gi.progbar.tics - 1); + double lo = gi.progbar.pos[gi.progbar.dim] - 0.5 * gi.progbar.length; + for (int i = 0; i < gi.progbar.tics; ++i) { + imgobjs[n] = DumpImage::CYLINDER; + imgparms[n][0] = gi.progbar.type1; + imgparms[n][1] = gi.progbar.pos[X]; + imgparms[n][2] = gi.progbar.pos[Y]; + imgparms[n][3] = gi.progbar.pos[Z]; + imgparms[n][4] = gi.progbar.pos[X]; + imgparms[n][5] = gi.progbar.pos[Y]; + imgparms[n][6] = gi.progbar.pos[Z]; + imgparms[n][7] = 1.1 * gi.progbar.diameter; + switch (gi.progbar.dim) { + case X: + imgparms[n][1] = lo + delta * i - 0.05 * delta; + imgparms[n][4] = lo + delta * i + 0.05 * delta; + break; + case Y: + imgparms[n][2] = lo + delta * i - 0.05 * delta; + imgparms[n][5] = lo + delta * i + 0.05 * delta; + break; + case Z: + imgparms[n][3] = lo + delta * i - 0.05 * delta; + imgparms[n][6] = lo + delta * i + 0.05 * delta; + break; + } + ++n; + } + if (gi.progbar.pstr) { + int ivar = input->variable->find(gi.progbar.pstr); + if (ivar < 0) + error->all(FLERR, Error::NOLASTLINE, "Variable name {} for fix graphics does not exist", + gi.progbar.pstr); + if (input->variable->equalstyle(ivar) == 0) + error->all(FLERR, Error::NOLASTLINE, + "Fix graphics variable {} is not equal-style variable", gi.progbar.pstr); + gi.progbar.pvar = ivar; + } + } + } + end_of_step(); +} + +/* ---------------------------------------------------------------------- */ + +void FixGraphics::end_of_step() +{ + // evaluate variable if necessary, wrap with clear/add + + if (varflag) modify->clearstep_compute(); + + int n = 0; + for (auto &gi : items) { + if (gi.style == SPHERE) { + if (gi.sphere.xstr) gi.sphere.pos[X] = input->variable->compute_equal(gi.sphere.xvar); + if (gi.sphere.ystr) gi.sphere.pos[Y] = input->variable->compute_equal(gi.sphere.yvar); + if (gi.sphere.zstr) gi.sphere.pos[Z] = input->variable->compute_equal(gi.sphere.zvar); + if (gi.sphere.dstr) gi.sphere.diameter = 2.0 * input->variable->compute_equal(gi.sphere.dvar); + imgparms[n][1] = gi.sphere.pos[X]; + imgparms[n][2] = gi.sphere.pos[Y]; + imgparms[n][3] = gi.sphere.pos[Z]; + imgparms[n][4] = gi.sphere.diameter; + ++n; + } else if (gi.style == CYLINDER) { + if (gi.cylinder.x1str) + gi.cylinder.pos1[X] = input->variable->compute_equal(gi.cylinder.x1var); + if (gi.cylinder.y1str) + gi.cylinder.pos1[Y] = input->variable->compute_equal(gi.cylinder.y1var); + if (gi.cylinder.z1str) + gi.cylinder.pos1[Z] = input->variable->compute_equal(gi.cylinder.z1var); + if (gi.cylinder.x2str) + gi.cylinder.pos2[X] = input->variable->compute_equal(gi.cylinder.x2var); + if (gi.cylinder.y2str) + gi.cylinder.pos2[Y] = input->variable->compute_equal(gi.cylinder.y2var); + if (gi.cylinder.z2str) + gi.cylinder.pos2[Z] = input->variable->compute_equal(gi.cylinder.z2var); + if (gi.cylinder.dstr) + gi.cylinder.diameter = 2.0 * input->variable->compute_equal(gi.cylinder.dvar); + imgparms[n][1] = gi.cylinder.pos1[X]; + imgparms[n][2] = gi.cylinder.pos1[Y]; + imgparms[n][3] = gi.cylinder.pos1[Z]; + imgparms[n][4] = gi.cylinder.pos2[X]; + imgparms[n][5] = gi.cylinder.pos2[Y]; + imgparms[n][6] = gi.cylinder.pos2[Z]; + imgparms[n][7] = gi.cylinder.diameter; + ++n; + } else if (gi.style == ARROW) { + if (gi.arrow.x1str) gi.arrow.tip[X] = input->variable->compute_equal(gi.arrow.x1var); + if (gi.arrow.y1str) gi.arrow.tip[Y] = input->variable->compute_equal(gi.arrow.y1var); + if (gi.arrow.z1str) gi.arrow.tip[Z] = input->variable->compute_equal(gi.arrow.z1var); + if (gi.arrow.x2str) gi.arrow.bot[X] = input->variable->compute_equal(gi.arrow.x2var); + if (gi.arrow.y2str) gi.arrow.bot[Y] = input->variable->compute_equal(gi.arrow.y2var); + if (gi.arrow.z2str) gi.arrow.bot[Z] = input->variable->compute_equal(gi.arrow.z2var); + if (gi.arrow.dstr) gi.arrow.diameter = 2.0 * input->variable->compute_equal(gi.arrow.dvar); + + double mid[3], vec[3]; + MathExtra::sub3(gi.arrow.bot, gi.arrow.tip, vec); + MathExtra::scaleadd3(gi.arrow.ratio, vec, gi.arrow.tip, mid); + imgparms[n][1] = gi.arrow.tip[X]; + imgparms[n][2] = gi.arrow.tip[Y]; + imgparms[n][3] = gi.arrow.tip[Z]; + imgparms[n][4] = mid[X]; + imgparms[n][5] = mid[Y]; + imgparms[n][6] = mid[Z]; + imgparms[n][7] = gi.arrow.diameter * (1.0 + 5.0 * gi.arrow.ratio); + ++n; + imgparms[n][1] = mid[X]; + imgparms[n][2] = mid[Y]; + imgparms[n][3] = mid[Z]; + imgparms[n][4] = gi.arrow.bot[X]; + imgparms[n][5] = gi.arrow.bot[Y]; + imgparms[n][6] = gi.arrow.bot[Z]; + imgparms[n][7] = gi.arrow.diameter; + ++n; + } else if (gi.style == PROGBAR) { + ++n; + if (gi.progbar.pstr) gi.progbar.progress = input->variable->compute_equal(gi.progbar.pvar); + // bracket into [0.0;1.0] rather than throwing an error for just a viz item + gi.progbar.progress = std::max(std::min(gi.progbar.progress, 1.0), 0.0); + switch (gi.progbar.dim) { + case X: + imgparms[n][1] = gi.progbar.pos[X] + (gi.progbar.progress - 0.5) * gi.progbar.length; + break; + case Y: + imgparms[n][2] = gi.progbar.pos[Y] + (gi.progbar.progress - 0.5) * gi.progbar.length; + break; + case Z: + imgparms[n][3] = gi.progbar.pos[Z] + (gi.progbar.progress - 0.5) * gi.progbar.length; + break; + default: + break; + } + ++n; + n += gi.progbar.tics; + } + } + if (varflag) modify->addstep_compute((update->ntimestep / nevery) * nevery + nevery); +} + +/* ---------------------------------------------------------------------- + provide graphics information to dump image +------------------------------------------------------------------------- */ + +int FixGraphics::image(int *&objs, double **&parms) +{ + objs = imgobjs; + parms = imgparms; + return numobjs; +} diff --git a/src/fix_graphics.h b/src/fix_graphics.h new file mode 100644 index 00000000000..0ec79ebdc21 --- /dev/null +++ b/src/fix_graphics.h @@ -0,0 +1,107 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS +// clang-format off +FixStyle(graphics,FixGraphics); +// clang-format on +#else + +#ifndef LMP_FIX_GRAPHICS_H +#define LMP_FIX_GRAPHICS_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixGraphics : public Fix { + public: + FixGraphics(class LAMMPS *, int, char **); + ~FixGraphics() override; + int setmask() override; + void init() override; + void end_of_step() override; + + int image(int *&, double **&) override; + + protected: + int varflag; + int numobjs; + int *imgobjs; + double **imgparms; + + struct SphereItem { + int style; + int type; + double pos[3]; + double diameter; + char *xstr, *ystr, *zstr, *dstr; + int xvar, yvar, zvar, dvar; + }; + + struct CylinderItem { + int style; + int type; + double pos1[3]; + double pos2[3]; + double diameter; + char *x1str, *y1str, *z1str, *x2str, *y2str, *z2str, *dstr; + int x1var, y1var, z1var, x2var, y2var, z2var, dvar; + }; + + struct ArrowItem { + int style; + int type; + double tip[3]; + double bot[3]; + double diameter; + double ratio; + char *x1str, *y1str, *z1str, *x2str, *y2str, *z2str, *dstr; + int x1var, y1var, z1var, x2var, y2var, z2var, dvar; + }; + + struct ProgbarItem { + int style; + int type1; + int type2; + int dim; + int tics; + double pos[3]; + double diameter; + double length; + double progress; + char *pstr; + int pvar; + }; + + union GraphicsItem { + GraphicsItem() = delete; + GraphicsItem(const SphereItem &s) : sphere(s) {} + GraphicsItem(const CylinderItem &c) : cylinder(c) {} + GraphicsItem(const ArrowItem &a) : arrow(a) {} + GraphicsItem(const ProgbarItem &p) : progbar(p) {} + + int style; + SphereItem sphere; + CylinderItem cylinder; + ArrowItem arrow; + ProgbarItem progbar; + }; + + std::vector items; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index f076d1a7492..0d98d43ba6f 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -19,10 +19,12 @@ #include "atom.h" #include "domain.h" +#include "dump_image.h" #include "error.h" #include "input.h" #include "lattice.h" #include "math_extra.h" +#include "memory.h" #include "modify.h" #include "respa.h" #include "update.h" @@ -41,7 +43,8 @@ enum { INSIDE, OUTSIDE }; FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr), pstr(nullptr), - rlostr(nullptr), rhistr(nullptr), lostr(nullptr), histr(nullptr) + rlostr(nullptr), rhistr(nullptr), lostr(nullptr), histr(nullptr), imgobjs(nullptr), + imgparms(nullptr) { if (narg < 4) utils::missing_cmd_args(FLERR, "fix indent", error); @@ -56,7 +59,7 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : ilevel_respa = 0; k = utils::numeric(FLERR, arg[3], false, lmp); - if (k < 0.0) error->all(FLERR, "Illegal fix indent force constant: {}", k); + if (k < 0.0) error->all(FLERR, 3, "Illegal fix indent force constant: {}", k); k3 = k / 3.0; // read geometry of indenter and optional args @@ -111,13 +114,45 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : pvalue *= zscale; } else - error->all(FLERR, "Unknown fix indent keyword: {}", istyle); + error->all(FLERR, "Unknown fix indent style: {}", istyle); varflag = 0; if (xstr || ystr || zstr || rstr || pstr || rlostr || rhistr || lostr || histr) varflag = 1; indenter_flag = 0; indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0; + + // set up indenter visualization + + if (istyle == SPHERE) { + // one sphere object to draw + memory->create(imgobjs, 1, "fix_indent:imgobjs"); + memory->create(imgparms, 1, 5, "fix_indent:imgparms"); + imgobjs[0] = DumpImage::SPHERE; + imgparms[0][0] = 1; // use color of first atom type + } else if (istyle == CYLINDER) { + // one cylinder object to draw + memory->create(imgobjs, 1, "fix_indent:imgobjs"); + memory->create(imgparms, 1, 8, "fix_indent:imgparms"); + imgobjs[0] = DumpImage::CYLINDER; + imgparms[0][0] = 1; // use color of first atom type + } else if (istyle == PLANE) { + if (domain->dimension == 2) { + // one cylinder object to draw in 2d + memory->create(imgobjs, 1, "fix_indent:imgobjs"); + memory->create(imgparms, 1, 8, "fix_indent:imgparms"); + imgobjs[0] = DumpImage::CYLINDER; + imgparms[0][0] = 1; // use color of first atom type + } else { + // two triangle objects to draw in 3d + memory->create(imgobjs, 2, "fix_indent:imgobjs"); + memory->create(imgparms, 2, 10, "fix_indent:imgparms"); + imgobjs[0] = DumpImage::TRIANGLE; + imgobjs[1] = DumpImage::TRIANGLE; + imgparms[0][0] = 1; // use color of first atom type by default + imgparms[1][0] = 1; // use color of first atom type by default + } + } } /* ---------------------------------------------------------------------- */ @@ -133,6 +168,9 @@ FixIndent::~FixIndent() delete[] rhistr; delete[] lostr; delete[] histr; + + memory->destroy(imgobjs); + memory->destroy(imgparms); } /* ---------------------------------------------------------------------- */ @@ -268,7 +306,7 @@ void FixIndent::post_force(int /*vflag*/) double radius = rstr ? input->variable->compute_equal(rvar) : rvalue; if (radius < 0.0) error->all(FLERR, "Illegal fix indent sphere radius: {}", radius); - for (int i = 0; i < nlocal; i++) + for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { delx = x[i][0] - ctr[0]; dely = x[i][1] - ctr[1]; @@ -294,6 +332,14 @@ void FixIndent::post_force(int /*vflag*/) indenter[2] -= fy; indenter[3] -= fz; } + } + + // store indenter object visualization parameters + + imgparms[0][1] = ctr[0]; + imgparms[0][2] = ctr[1]; + imgparms[0][3] = ctr[2]; + imgparms[0][4] = 2.0 * radius; // cylindrical indenter @@ -309,7 +355,7 @@ void FixIndent::post_force(int /*vflag*/) double radius{rstr ? input->variable->compute_equal(rvar) : rvalue}; if (radius < 0.0) error->all(FLERR, "Illegal fix indent cylinder radius: {}", radius); - for (int i = 0; i < nlocal; i++) + for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { double del[3] = {x[i][0] - ctr[0], x[i][1] - ctr[1], x[i][2] - ctr[2]}; del[cdim] = 0; @@ -334,6 +380,18 @@ void FixIndent::post_force(int /*vflag*/) indenter[2] -= fy; indenter[3] -= fz; } + } + + // store indenter object visualization parameters: positions of cylinder edges and diameter + + imgparms[0][1] = ctr[0]; + imgparms[0][2] = ctr[1]; + imgparms[0][3] = ctr[2]; + ctr[cdim] = domain->boxhi[cdim]; + imgparms[0][4] = ctr[0]; + imgparms[0][5] = ctr[1]; + imgparms[0][6] = ctr[2]; + imgparms[0][7] = 2.0 * radius; // conical indenter @@ -403,13 +461,13 @@ void FixIndent::post_force(int /*vflag*/) // planar indenter - } else { + } else { // if (istyle == PLANE) // plane = current plane position double plane{pstr ? input->variable->compute_equal(pvar) : pvalue}; - for (int i = 0; i < nlocal; i++) + for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { dr = planeside * (plane - x[i][cdim]); if (dr >= 0.0) continue; @@ -418,6 +476,98 @@ void FixIndent::post_force(int /*vflag*/) indenter[0] -= k3 * dr * dr * dr; indenter[cdim + 1] -= fmag; } + } + + // store indenter object visualization parameters + + if (domain->dimension == 2) { + switch (cdim) { + case 0: + imgparms[0][1] = planeside * plane; + imgparms[0][2] = domain->boxlo[1]; + imgparms[0][3] = 0.0; + imgparms[0][4] = planeside * plane; + imgparms[0][5] = domain->boxhi[1]; + imgparms[0][6] = 0.0; + imgparms[0][7] = 0.0; // no simple guess for diameter. need to use fflag2 to adjust + break; + case 1: + imgparms[0][1] = domain->boxlo[0]; + imgparms[0][2] = planeside * plane; + imgparms[0][3] = 0.0; + imgparms[0][4] = domain->boxhi[0]; + imgparms[0][5] = planeside * plane; + imgparms[0][6] = 0.0; + imgparms[0][7] = 0.0; // no simple guess for diameter. need to use fflag2 to adjust + break; + case 2:; // no planar indenter allowed in z-direction for 2d systems + break; + } + } else { + // two triangles + switch (cdim) { + case 0: + imgparms[0][1] = planeside * plane; + imgparms[0][2] = domain->boxlo[1]; + imgparms[0][3] = domain->boxlo[2]; + imgparms[0][4] = planeside * plane; + imgparms[0][5] = domain->boxhi[1]; + imgparms[0][6] = domain->boxlo[2]; + imgparms[0][7] = planeside * plane; + imgparms[0][8] = domain->boxlo[1]; + imgparms[0][9] = domain->boxhi[2]; + imgparms[1][1] = planeside * plane; + imgparms[1][2] = domain->boxhi[1]; + imgparms[1][3] = domain->boxhi[2]; + imgparms[1][4] = planeside * plane; + imgparms[1][5] = domain->boxlo[1]; + imgparms[1][6] = domain->boxhi[2]; + imgparms[1][7] = planeside * plane; + imgparms[1][8] = domain->boxhi[1]; + imgparms[1][9] = domain->boxlo[2]; + break; + case 1: + imgparms[0][1] = domain->boxlo[0]; + imgparms[0][2] = planeside * plane; + imgparms[0][3] = domain->boxlo[2]; + imgparms[0][4] = domain->boxhi[0]; + imgparms[0][5] = planeside * plane; + imgparms[0][6] = domain->boxlo[2]; + imgparms[0][7] = domain->boxlo[0]; + imgparms[0][8] = planeside * plane; + imgparms[0][9] = domain->boxhi[2]; + imgparms[1][1] = domain->boxhi[0]; + imgparms[1][2] = planeside * plane; + imgparms[1][3] = domain->boxhi[2]; + imgparms[1][4] = domain->boxlo[0]; + imgparms[1][5] = planeside * plane; + imgparms[1][6] = domain->boxhi[2]; + imgparms[1][7] = domain->boxhi[0]; + imgparms[1][8] = planeside * plane; + imgparms[1][9] = domain->boxlo[2]; + break; + case 2: + imgparms[0][1] = domain->boxlo[0]; + imgparms[0][2] = domain->boxlo[1]; + imgparms[0][3] = planeside * plane; + imgparms[0][4] = domain->boxhi[0]; + imgparms[0][5] = domain->boxlo[1]; + imgparms[0][6] = planeside * plane; + imgparms[0][7] = domain->boxlo[0]; + imgparms[0][8] = domain->boxhi[1]; + imgparms[0][9] = planeside * plane; + imgparms[1][1] = domain->boxhi[0]; + imgparms[1][2] = domain->boxhi[1]; + imgparms[1][3] = planeside * plane; + imgparms[1][4] = domain->boxlo[0]; + imgparms[1][5] = domain->boxhi[1]; + imgparms[1][6] = planeside * plane; + imgparms[1][7] = domain->boxhi[0]; + imgparms[1][8] = domain->boxlo[1]; + imgparms[1][9] = planeside * plane; + break; + } + } } if (varflag) modify->addstep_compute(update->ntimestep + 1); @@ -853,3 +1003,24 @@ double FixIndent::closest(double *x, double *near, double *nearest, double dsq) nearest[2] = near[2]; return rsq; } + +/* ---------------------------------------------------------------------- + provide graphics information to dump image to render indenter + data has been copied to dedicated storage during fix indent execution +------------------------------------------------------------------------- */ +int FixIndent::image(int *&objs, double **&parms) +{ + objs = imgobjs; + parms = imgparms; + if (istyle == SPHERE) + return 1; + else if (istyle == CYLINDER) + return 1; + else if (istyle == PLANE) + if (domain->dimension == 2) + return 1; + else + return 2; + else + return 0; +} diff --git a/src/fix_indent.h b/src/fix_indent.h index 833f576b467..01b54657673 100644 --- a/src/fix_indent.h +++ b/src/fix_indent.h @@ -38,6 +38,8 @@ class FixIndent : public Fix { double compute_scalar() override; double compute_vector(int) override; + int image(int *&, double **&) override; + private: int istyle, scaleflag, side; double k, k3; @@ -53,6 +55,11 @@ class FixIndent : public Fix { int rlovar, rhivar, lovar, hivar; double rlovalue, rhivalue, lovalue, hivalue; + // arrays for dump image rendering + + int *imgobjs; + double **imgparms; + // methods for argument parsing int geometry(int, char **); diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp index 303994c5736..93ebc0c0fc4 100644 --- a/src/fix_wall.cpp +++ b/src/fix_wall.cpp @@ -14,9 +14,11 @@ #include "fix_wall.h" #include "domain.h" +#include "dump_image.h" #include "error.h" #include "input.h" #include "lattice.h" +#include "memory.h" #include "modify.h" #include "respa.h" #include "update.h" @@ -29,6 +31,107 @@ using namespace FixConst; static const char *wallpos[] = {"xlo", "xhi", "ylo", "yhi", "zlo", "zhi"}; +// record wall info for dump image + +void FixWall::update_image_plane(int m, int which, double coord, double **imgparms, Domain *domain) +{ + if (domain->dimension == 2) { + // one cylinder for 2d. diameter is zero and can be set with fparam2 + switch (which) { + case XLO: // fallthrough + case XHI: + imgparms[m][1] = coord; + imgparms[m][2] = domain->boxlo[1]; + imgparms[m][3] = 0.0; + imgparms[m][4] = coord; + imgparms[m][5] = domain->boxhi[1]; + imgparms[m][6] = 0.0; + imgparms[m][7] = 0.0; + break; + case YLO: // fallthrough + case YHI: + imgparms[m][1] = domain->boxlo[0]; + imgparms[m][2] = coord; + imgparms[m][3] = 0.0; + imgparms[m][4] = domain->boxhi[0]; + imgparms[m][5] = coord; + imgparms[m][6] = 0.0; + imgparms[m][7] = 0.0; + break; + case ZLO: // fallthrough + case ZHI:; // no wall in z-direction allowed for 2d systems + break; + } + } else { + // two triangles for 3d + switch (which) { + case XLO: // fallthrough + case XHI: + imgparms[2 * m][1] = coord; + imgparms[2 * m][2] = domain->boxlo[1]; + imgparms[2 * m][3] = domain->boxlo[2]; + imgparms[2 * m][4] = coord; + imgparms[2 * m][5] = domain->boxhi[1]; + imgparms[2 * m][6] = domain->boxlo[2]; + imgparms[2 * m][7] = coord; + imgparms[2 * m][8] = domain->boxlo[1]; + imgparms[2 * m][9] = domain->boxhi[2]; + imgparms[2 * m + 1][1] = coord; + imgparms[2 * m + 1][2] = domain->boxhi[1]; + imgparms[2 * m + 1][3] = domain->boxhi[2]; + imgparms[2 * m + 1][4] = coord; + imgparms[2 * m + 1][5] = domain->boxlo[1]; + imgparms[2 * m + 1][6] = domain->boxhi[2]; + imgparms[2 * m + 1][7] = coord; + imgparms[2 * m + 1][8] = domain->boxhi[1]; + imgparms[2 * m + 1][9] = domain->boxlo[2]; + break; + case YLO: // fallthrough + case YHI: + imgparms[2 * m][1] = domain->boxlo[0]; + imgparms[2 * m][2] = coord; + imgparms[2 * m][3] = domain->boxlo[2]; + imgparms[2 * m][4] = domain->boxhi[0]; + imgparms[2 * m][5] = coord; + imgparms[2 * m][6] = domain->boxlo[2]; + imgparms[2 * m][7] = domain->boxlo[0]; + imgparms[2 * m][8] = coord; + imgparms[2 * m][9] = domain->boxhi[2]; + imgparms[2 * m + 1][1] = domain->boxhi[0]; + imgparms[2 * m + 1][2] = coord; + imgparms[2 * m + 1][3] = domain->boxhi[2]; + imgparms[2 * m + 1][4] = domain->boxlo[0]; + imgparms[2 * m + 1][5] = coord; + imgparms[2 * m + 1][6] = domain->boxhi[2]; + imgparms[2 * m + 1][7] = domain->boxhi[0]; + imgparms[2 * m + 1][8] = coord; + imgparms[2 * m + 1][9] = domain->boxlo[2]; + break; + case ZLO: // fallthrough + case ZHI: + imgparms[2 * m][1] = domain->boxlo[0]; + imgparms[2 * m][2] = domain->boxlo[1]; + imgparms[2 * m][3] = coord; + imgparms[2 * m][4] = domain->boxhi[0]; + imgparms[2 * m][5] = domain->boxlo[1]; + imgparms[2 * m][6] = coord; + imgparms[2 * m][7] = domain->boxlo[0]; + imgparms[2 * m][8] = domain->boxhi[1]; + imgparms[2 * m][9] = coord; + imgparms[2 * m + 1][1] = domain->boxhi[0]; + imgparms[2 * m + 1][2] = domain->boxhi[1]; + imgparms[2 * m + 1][3] = coord; + imgparms[2 * m + 1][4] = domain->boxlo[0]; + imgparms[2 * m + 1][5] = domain->boxhi[1]; + imgparms[2 * m + 1][6] = coord; + imgparms[2 * m + 1][7] = domain->boxhi[0]; + imgparms[2 * m + 1][8] = domain->boxlo[1]; + imgparms[2 * m + 1][9] = coord; + break; + } + } +} + /* ---------------------------------------------------------------------- */ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), nwall(0) @@ -232,6 +335,27 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), nwall eflag = 0; for (int m = 0; m <= nwall; m++) ewall[m] = 0.0; + + // for rendering walls with dump image. + if (domain->dimension == 2) { + // one cylinder object per wall to draw in 2d + memory->create(imgobjs, nwall, "fix_wall:imgobjs"); + memory->create(imgparms, nwall, 8, "fix_wall:imgparms"); + for (int m = 0; m < nwall; ++m) { + imgobjs[m] = DumpImage::CYLINDER; + imgparms[m][0] = 1; // use color of first atom type by default + } + } else { + // two triangle objects per wall to draw in 3d + memory->create(imgobjs, 2 * nwall, "fix_wall:imgobjs"); + memory->create(imgparms, 2 * nwall, 10, "fix_wall:imgparms"); + for (int m = 0; m < nwall; ++m) { + imgobjs[2 * m] = DumpImage::TRIANGLE; + imgobjs[2 * m + 1] = DumpImage::TRIANGLE; + imgparms[2 * m][0] = 1; // use color of first atom type by default + imgparms[2 * m + 1][0] = 1; // use color of first atom type by default + } + } } /* ---------------------------------------------------------------------- */ @@ -248,6 +372,9 @@ FixWall::~FixWall() delete[] fstr[m]; delete[] kstr[m]; } + + memory->destroy(imgobjs); + memory->destroy(imgparms); } /* ---------------------------------------------------------------------- */ @@ -325,7 +452,7 @@ void FixWall::min_setup(int vflag) /* ---------------------------------------------------------------------- only called if fldflag set, in place of post_force -------------------------------------------------------------------------- */ + ------------------------------------------------------------------------- */ void FixWall::pre_force(int vflag) { @@ -377,8 +504,9 @@ void FixWall::post_force(int vflag) } wall_particle(m, wallwhich[m], coord); - } + update_image_plane(m, wallwhich[m], coord, imgparms, domain); + } if (varflag) modify->addstep_compute(update->ntimestep + 1); } @@ -425,3 +553,19 @@ double FixWall::compute_vector(int n) } return ewall_all[n + 1]; } + +/* ---------------------------------------------------------------------- + provide graphics information to dump image to render wall as plane + data has been copied to dedicated storage during fix indent execution +------------------------------------------------------------------------- */ + +int FixWall::image(int *&objs, double **&parms) +{ + objs = imgobjs; + parms = imgparms; + if (domain->dimension == 2) { + return nwall; + } else { + return 2 * nwall; + } +} diff --git a/src/fix_wall.h b/src/fix_wall.h index 871f4acf495..682d9785e3e 100644 --- a/src/fix_wall.h +++ b/src/fix_wall.h @@ -43,8 +43,11 @@ class FixWall : public Fix { double compute_scalar() override; double compute_vector(int) override; + int image(int *&, double **&) override; + virtual void precompute(int) = 0; virtual void wall_particle(int, int, double) = 0; + static void update_image_plane(int, int, double, double **, class Domain *); protected: double epsilon[6], sigma[6], alpha[6], cutoff[6]; @@ -57,6 +60,8 @@ class FixWall : public Fix { int eflag; // per-wall flag for energy summation int ilevel_respa; int fldflag; + int *imgobjs; + double **imgparms; }; } // namespace LAMMPS_NS diff --git a/src/fix_wall_reflect.cpp b/src/fix_wall_reflect.cpp index 0169644e4ac..43b9bd51f44 100644 --- a/src/fix_wall_reflect.cpp +++ b/src/fix_wall_reflect.cpp @@ -17,9 +17,12 @@ #include "atom.h" #include "comm.h" #include "domain.h" +#include "dump_image.h" #include "error.h" +#include "fix_wall.h" #include "input.h" #include "lattice.h" +#include "memory.h" #include "modify.h" #include "update.h" #include "variable.h" @@ -32,7 +35,7 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), nwall(0), varflag(0) + Fix(lmp, narg, arg), nwall(0), varflag(0), imgobjs(nullptr), imgparms(nullptr) { if (narg < 4) utils::missing_cmd_args(FLERR, "fix wall/reflect", error); @@ -55,12 +58,12 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : if (iarg+2 > narg) error->all(FLERR, "Illegal fix wall/reflect {} command: missing argument(s)", arg[iarg]); int newwall; - if (strcmp(arg[iarg],"xlo") == 0) newwall = XLO; - else if (strcmp(arg[iarg],"xhi") == 0) newwall = XHI; - else if (strcmp(arg[iarg],"ylo") == 0) newwall = YLO; - else if (strcmp(arg[iarg],"yhi") == 0) newwall = YHI; - else if (strcmp(arg[iarg],"zlo") == 0) newwall = ZLO; - else if (strcmp(arg[iarg],"zhi") == 0) newwall = ZHI; + if (strcmp(arg[iarg],"xlo") == 0) newwall = FixWall::XLO; + else if (strcmp(arg[iarg],"xhi") == 0) newwall = FixWall::XHI; + else if (strcmp(arg[iarg],"ylo") == 0) newwall = FixWall::YLO; + else if (strcmp(arg[iarg],"yhi") == 0) newwall = FixWall::YHI; + else if (strcmp(arg[iarg],"zlo") == 0) newwall = FixWall::ZLO; + else if (strcmp(arg[iarg],"zhi") == 0) newwall = FixWall::ZHI; for (int m = 0; (m < nwall) && (m < 6); m++) if (newwall == wallwhich[m]) @@ -99,16 +102,16 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : if (nwall == 0) utils::missing_cmd_args(FLERR, "fix wall/reflect", error); for (int m = 0; m < nwall; m++) { - if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic) + if ((wallwhich[m] == FixWall::XLO || wallwhich[m] == FixWall::XHI) && domain->xperiodic) error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension x"); - if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic) + if ((wallwhich[m] == FixWall::YLO || wallwhich[m] == FixWall::YHI) && domain->yperiodic) error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension y"); - if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic) + if ((wallwhich[m] == FixWall::ZLO || wallwhich[m] == FixWall::ZHI) && domain->zperiodic) error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension z"); } for (int m = 0; m < nwall; m++) - if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2) + if ((wallwhich[m] == FixWall::ZLO || wallwhich[m] == FixWall::ZHI) && domain->dimension == 2) error->all(FLERR, "Cannot use fix wall/reflect zlo/zhi for a 2d simulation"); @@ -128,8 +131,8 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : for (int m = 0; m < nwall; m++) { if (wallstyle[m] != CONSTANT) continue; - if (wallwhich[m] < YLO) coord0[m] *= xscale; - else if (wallwhich[m] < ZLO) coord0[m] *= yscale; + if (wallwhich[m] < FixWall::YLO) coord0[m] *= xscale; + else if (wallwhich[m] < FixWall::ZLO) coord0[m] *= yscale; else coord0[m] *= zscale; } } @@ -139,6 +142,27 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : varflag = 0; for (int m = 0; m < nwall; m++) if (wallstyle[m] == VARIABLE) varflag = 1; + + // for rendering walls with dump image. + if (domain->dimension == 2) { + // one cylinder object per wall to draw in 2d + memory->create(imgobjs, nwall, "fix_wall_reflect:imgobjs"); + memory->create(imgparms, nwall, 8, "fix_wall_reflect:imgparms"); + for (int m = 0; m < nwall; ++m) { + imgobjs[m] = DumpImage::CYLINDER; + imgparms[m][0] = 1; // use color of first atom type by default + } + } else { + // two triangle objects per wall to draw in 3d + memory->create(imgobjs, 2 * nwall, "fix_wall_reflect:imgobjs"); + memory->create(imgparms, 2 * nwall, 10, "fix_wall_reflect:imgparms"); + for (int m = 0; m < nwall; ++m) { + imgobjs[2 * m] = DumpImage::TRIANGLE; + imgobjs[2 * m + 1] = DumpImage::TRIANGLE; + imgparms[2 * m][0] = 1; // use color of first atom type by default + imgparms[2 * m + 1][0] = 1; // use color of first atom type by default + } + } } /* ---------------------------------------------------------------------- */ @@ -149,6 +173,9 @@ FixWallReflect::~FixWallReflect() for (int m = 0; m < nwall; m++) if (wallstyle[m] == VARIABLE) delete [] varstr[m]; + + memory->destroy(imgobjs); + memory->destroy(imgparms); } /* ---------------------------------------------------------------------- */ @@ -196,12 +223,15 @@ void FixWallReflect::post_integrate() for (int m = 0; m < nwall; m++) { if (wallstyle[m] == VARIABLE) { coord = input->variable->compute_equal(varindex[m]); - if (wallwhich[m] < YLO) coord *= xscale; - else if (wallwhich[m] < ZLO) coord *= yscale; + if (wallwhich[m] < FixWall::YLO) coord *= xscale; + else if (wallwhich[m] < FixWall::ZLO) coord *= yscale; else coord *= zscale; } else coord = coord0[m]; wall_particle(m,wallwhich[m],coord); + + // record wall graphics objects for dump image + FixWall::update_image_plane(m, wallwhich[m], coord, imgparms, domain); } if (varflag) modify->addstep_compute(update->ntimestep + 1); @@ -239,3 +269,19 @@ void FixWallReflect::wall_particle(int /* m */, int which, double coord) } } } + +/* ---------------------------------------------------------------------- + provide graphics information to dump image to render wall as plane + data has been copied to dedicated storage during fix indent execution +------------------------------------------------------------------------- */ + +int FixWallReflect::image(int *&objs, double **&parms) +{ + objs = imgobjs; + parms = imgparms; + if (domain->dimension == 2) { + return nwall; + } else { + return 2 * nwall; + } +} diff --git a/src/fix_wall_reflect.h b/src/fix_wall_reflect.h index 05acfe94cf1..8f8e8efd07d 100644 --- a/src/fix_wall_reflect.h +++ b/src/fix_wall_reflect.h @@ -26,7 +26,6 @@ namespace LAMMPS_NS { class FixWallReflect : public Fix { public: - enum { XLO = 0, XHI = 1, YLO = 2, YHI = 3, ZLO = 4, ZHI = 5 }; enum { NONE = 0, EDGE, CONSTANT, VARIABLE }; FixWallReflect(class LAMMPS *, int, char **); @@ -35,6 +34,8 @@ class FixWallReflect : public Fix { void init() override; void post_integrate() override; + int image(int *&, double **&) override; + protected: int nwall; int wallwhich[6], wallstyle[6]; @@ -43,6 +44,8 @@ class FixWallReflect : public Fix { int varindex[6]; int varflag; double xscale, yscale, zscale; + int *imgobjs; + double **imgparms; virtual void wall_particle(int m, int which, double coord); }; diff --git a/src/image.cpp b/src/image.cpp index c806268617f..c8ad2d351fb 100644 --- a/src/image.cpp +++ b/src/image.cpp @@ -44,20 +44,120 @@ using MathConst::DEG2RAD; using MathConst::MY_PI; using MathConst::MY_PI4; -static constexpr int NCOLORS = 140; -static constexpr int NELEMENTS = 109; -static constexpr double EPSILON = 1.0e-6; - -enum{NUMERIC,MINVALUE,MAXVALUE}; -enum{CONTINUOUS,DISCRETE,SEQUENTIAL}; -enum{ABSOLUTE,FRACTIONAL}; -enum{NO,YES}; +// clang-format on +namespace { +constexpr int NCOLORS = 140; +constexpr int NELEMENTS = 109; +constexpr double EPSILON = 1.0e-6; + +enum { NUMERIC, MINVALUE, MAXVALUE }; +enum { CONTINUOUS, DISCRETE, SEQUENTIAL }; +enum { ABSOLUTE, FRACTIONAL }; +enum { NO, YES }; + +//////////////////////////////////////////////////////////////////////// +// the following regular Bayer threshold matrix can be created for any +// power of 2 ranks with the following python code. +// See https://en.wikipedia.org/wiki/Ordered_dithering +// +// import numpy as np +// +// def bayer_matrix(n: int) -> np.ndarray: +// """Generate an n x n Bayer (ordered dither) matrix for n a power of 2.""" +// if n == 1: +// return np.array([[0]], dtype=int) +// +// m = n // 2 +// B = bayer_matrix(m) +// +// # Standard Bayer recursion +// return np.block([ +// [4 * B + 0, 4 * B + 2], +// [4 * B + 3, 4 * B + 1] +// ]) +// +// # set rank +// n=16 +// +// # normalized matrix from recursion +// matrix = (bayer_matrix(n) + 0.5) / (n * n) +// +// print("constexpr int TRANK = %d;" % n) +// print("constexpr double transthresh[TRANK][TRANK] = {") +// nx, ny = matrix.shape +// for iy in range(0,ny): +// print("{") +// for ix in range(0,nx): +// if ix < nx -1: +// print(f"{matrix[ix][iy]:12.9f},") +// else: +// print(f"{matrix[ix][iy]:12.9f}") +// if iy < ny - 1: +// print("},") +// else: +// print("}") +// +// print("};") +//////////////////////////////////////////////////////////////////////// +constexpr int TRANK = 16; +constexpr double transthresh[TRANK][TRANK] = { + {0.001953125, 0.751953125, 0.189453125, 0.939453125, 0.048828125, 0.798828125, 0.236328125, + 0.986328125, 0.013671875, 0.763671875, 0.201171875, 0.951171875, 0.060546875, 0.810546875, + 0.248046875, 0.998046875}, + {0.501953125, 0.251953125, 0.689453125, 0.439453125, 0.548828125, 0.298828125, 0.736328125, + 0.486328125, 0.513671875, 0.263671875, 0.701171875, 0.451171875, 0.560546875, 0.310546875, + 0.748046875, 0.498046875}, + {0.126953125, 0.876953125, 0.064453125, 0.814453125, 0.173828125, 0.923828125, 0.111328125, + 0.861328125, 0.138671875, 0.888671875, 0.076171875, 0.826171875, 0.185546875, 0.935546875, + 0.123046875, 0.873046875}, + {0.626953125, 0.376953125, 0.564453125, 0.314453125, 0.673828125, 0.423828125, 0.611328125, + 0.361328125, 0.638671875, 0.388671875, 0.576171875, 0.326171875, 0.685546875, 0.435546875, + 0.623046875, 0.373046875}, + {0.033203125, 0.783203125, 0.220703125, 0.970703125, 0.017578125, 0.767578125, 0.205078125, + 0.955078125, 0.044921875, 0.794921875, 0.232421875, 0.982421875, 0.029296875, 0.779296875, + 0.216796875, 0.966796875}, + {0.533203125, 0.283203125, 0.720703125, 0.470703125, 0.517578125, 0.267578125, 0.705078125, + 0.455078125, 0.544921875, 0.294921875, 0.732421875, 0.482421875, 0.529296875, 0.279296875, + 0.716796875, 0.466796875}, + {0.158203125, 0.908203125, 0.095703125, 0.845703125, 0.142578125, 0.892578125, 0.080078125, + 0.830078125, 0.169921875, 0.919921875, 0.107421875, 0.857421875, 0.154296875, 0.904296875, + 0.091796875, 0.841796875}, + {0.658203125, 0.408203125, 0.595703125, 0.345703125, 0.642578125, 0.392578125, 0.580078125, + 0.330078125, 0.669921875, 0.419921875, 0.607421875, 0.357421875, 0.654296875, 0.404296875, + 0.591796875, 0.341796875}, + {0.009765625, 0.759765625, 0.197265625, 0.947265625, 0.056640625, 0.806640625, 0.244140625, + 0.994140625, 0.005859375, 0.755859375, 0.193359375, 0.943359375, 0.052734375, 0.802734375, + 0.240234375, 0.990234375}, + {0.509765625, 0.259765625, 0.697265625, 0.447265625, 0.556640625, 0.306640625, 0.744140625, + 0.494140625, 0.505859375, 0.255859375, 0.693359375, 0.443359375, 0.552734375, 0.302734375, + 0.740234375, 0.490234375}, + {0.134765625, 0.884765625, 0.072265625, 0.822265625, 0.181640625, 0.931640625, 0.119140625, + 0.869140625, 0.130859375, 0.880859375, 0.068359375, 0.818359375, 0.177734375, 0.927734375, + 0.115234375, 0.865234375}, + {0.634765625, 0.384765625, 0.572265625, 0.322265625, 0.681640625, 0.431640625, 0.619140625, + 0.369140625, 0.630859375, 0.380859375, 0.568359375, 0.318359375, 0.677734375, 0.427734375, + 0.615234375, 0.365234375}, + {0.041015625, 0.791015625, 0.228515625, 0.978515625, 0.025390625, 0.775390625, 0.212890625, + 0.962890625, 0.037109375, 0.787109375, 0.224609375, 0.974609375, 0.021484375, 0.771484375, + 0.208984375, 0.958984375}, + {0.541015625, 0.291015625, 0.728515625, 0.478515625, 0.525390625, 0.275390625, 0.712890625, + 0.462890625, 0.537109375, 0.287109375, 0.724609375, 0.474609375, 0.521484375, 0.271484375, + 0.708984375, 0.458984375}, + {0.166015625, 0.916015625, 0.103515625, 0.853515625, 0.150390625, 0.900390625, 0.087890625, + 0.837890625, 0.162109375, 0.912109375, 0.099609375, 0.849609375, 0.146484375, 0.896484375, + 0.083984375, 0.833984375}, + {0.666015625, 0.416015625, 0.603515625, 0.353515625, 0.650390625, 0.400390625, 0.587890625, + 0.337890625, 0.662109375, 0.412109375, 0.599609375, 0.349609375, 0.646484375, 0.396484375, + 0.583984375, 0.333984375}}; +} // namespace +// clang-format off /* ---------------------------------------------------------------------- */ Image::Image(LAMMPS *lmp, int nmap_caller) : - Pointers(lmp), depthBuffer(nullptr), surfaceBuffer(nullptr), depthcopy(nullptr), - surfacecopy(nullptr), imageBuffer(nullptr), rgbcopy(nullptr), writeBuffer(nullptr) + Pointers(lmp), maps(nullptr), depthBuffer(nullptr), surfaceBuffer(nullptr), depthcopy(nullptr), + surfacecopy(nullptr), imageBuffer(nullptr), rgbcopy(nullptr), writeBuffer(nullptr), + recvcounts(nullptr), displs(nullptr), username(nullptr), userrgb(nullptr), random(nullptr) { MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); @@ -80,9 +180,6 @@ Image::Image(LAMMPS *lmp, int nmap_caller) : // colors ncolors = 0; - username = nullptr; - userrgb = nullptr; - boxcolor = color2rgb("yellow"); background[0] = background[1] = background[2] = 0; @@ -117,13 +214,6 @@ Image::Image(LAMMPS *lmp, int nmap_caller) : backLightColor[0] = 0.9; backLightColor[1] = 0.9; backLightColor[2] = 0.9; - - random = nullptr; - - // MPI_Gatherv vectors - - recvcounts = nullptr; - displs = nullptr; } /* ---------------------------------------------------------------------- */ @@ -131,7 +221,7 @@ Image::Image(LAMMPS *lmp, int nmap_caller) : Image::~Image() { for (int i = 0; i < nmap; i++) delete maps[i]; - delete [] maps; + delete[] maps; for (int i = 0; i < ncolors; i++) delete [] username[i]; memory->sfree(username); @@ -144,7 +234,7 @@ Image::~Image() memory->destroy(surfacecopy); memory->destroy(rgbcopy); - if (random) delete random; + delete random; memory->destroy(recvcounts); memory->destroy(displs); @@ -416,20 +506,20 @@ void Image::merge() draw simulation bounding box as 12 cylinders ------------------------------------------------------------------------- */ -void Image::draw_box(double (*corners)[3], double diameter) +void Image::draw_box(double (*corners)[3], double diameter, double opacity) { - draw_cylinder(corners[0],corners[1],boxcolor,diameter,3); - draw_cylinder(corners[2],corners[3],boxcolor,diameter,3); - draw_cylinder(corners[0],corners[2],boxcolor,diameter,3); - draw_cylinder(corners[1],corners[3],boxcolor,diameter,3); - draw_cylinder(corners[0],corners[4],boxcolor,diameter,3); - draw_cylinder(corners[1],corners[5],boxcolor,diameter,3); - draw_cylinder(corners[2],corners[6],boxcolor,diameter,3); - draw_cylinder(corners[3],corners[7],boxcolor,diameter,3); - draw_cylinder(corners[4],corners[5],boxcolor,diameter,3); - draw_cylinder(corners[6],corners[7],boxcolor,diameter,3); - draw_cylinder(corners[4],corners[6],boxcolor,diameter,3); - draw_cylinder(corners[5],corners[7],boxcolor,diameter,3); + draw_cylinder(corners[0],corners[1],boxcolor,diameter,3,opacity); + draw_cylinder(corners[2],corners[3],boxcolor,diameter,3,opacity); + draw_cylinder(corners[0],corners[2],boxcolor,diameter,3,opacity); + draw_cylinder(corners[1],corners[3],boxcolor,diameter,3,opacity); + draw_cylinder(corners[0],corners[4],boxcolor,diameter,3,opacity); + draw_cylinder(corners[1],corners[5],boxcolor,diameter,3,opacity); + draw_cylinder(corners[2],corners[6],boxcolor,diameter,3,opacity); + draw_cylinder(corners[3],corners[7],boxcolor,diameter,3,opacity); + draw_cylinder(corners[4],corners[5],boxcolor,diameter,3,opacity); + draw_cylinder(corners[6],corners[7],boxcolor,diameter,3,opacity); + draw_cylinder(corners[4],corners[6],boxcolor,diameter,3,opacity); + draw_cylinder(corners[5],corners[7],boxcolor,diameter,3,opacity); } /* ---------------------------------------------------------------------- @@ -437,11 +527,11 @@ void Image::draw_box(double (*corners)[3], double diameter) axes = 4 end points ------------------------------------------------------------------------- */ -void Image::draw_axes(double (*axes)[3], double diameter) +void Image::draw_axes(double (*axes)[3], double diameter, double opacity) { - draw_cylinder(axes[0],axes[1],color2rgb("red"),diameter,3); - draw_cylinder(axes[0],axes[2],color2rgb("green"),diameter,3); - draw_cylinder(axes[0],axes[3],color2rgb("blue"),diameter,3); + draw_cylinder(axes[0],axes[1],color2rgb("red"),diameter,3,opacity); + draw_cylinder(axes[0],axes[2],color2rgb("green"),diameter,3,opacity); + draw_cylinder(axes[0],axes[3],color2rgb("blue"),diameter,3,opacity); } /* ---------------------------------------------------------------------- @@ -449,7 +539,7 @@ void Image::draw_axes(double (*axes)[3], double diameter) render pixel by pixel onto image plane with depth buffering ------------------------------------------------------------------------- */ -void Image::draw_sphere(const double *x, const double *surfaceColor, double diameter) +void Image::draw_sphere(const double *x, const double *surfaceColor, double diameter, double opacity) { double xlocal[3]; @@ -483,8 +573,9 @@ void Image::draw_sphere(const double *x, const double *surfaceColor, double diam for (int iy = yc - pixelRadius; iy <= yc + pixelRadius; iy++) { for (int ix = xc - pixelRadius; ix <= xc + pixelRadius; ix++) { if (iy < 0 || iy >= height || ix < 0 || ix >= width) continue; - double surface[3]; + if (((opacity < 1.0) && (transthresh[ix % TRANK][iy % TRANK] > opacity)) || (opacity <= 0.0)) continue; + double surface[3]; surface[1] = ((iy - yc) - height_error) * pixelWidth; surface[0] = ((ix - xc) - width_error) * pixelWidth; double projRad = surface[0]*surface[0] + surface[1]*surface[1]; @@ -509,7 +600,7 @@ void Image::draw_sphere(const double *x, const double *surfaceColor, double diam render pixel by pixel onto image plane with depth buffering ------------------------------------------------------------------------- */ -void Image::draw_cube(const double *x, const double *surfaceColor, double diameter) +void Image::draw_cube(const double *x, const double *surfaceColor, double diameter, double opacity) { double xlocal[3],surface[3]; double normal[3] = {0.0, 0.0, 1.0}; @@ -548,6 +639,7 @@ void Image::draw_cube(const double *x, const double *surfaceColor, double diamet for (int iy = yc - pixelHalfWidth; iy <= yc + pixelHalfWidth; iy ++) { for (int ix = xc - pixelHalfWidth; ix <= xc + pixelHalfWidth; ix ++) { if (iy < 0 || iy >= height || ix < 0 || ix >= width) continue; + if (((opacity < 1.0) && (transthresh[ix % TRANK][iy % TRANK] > opacity)) || (opacity <= 0.0)) continue; double sy = ((iy - yc) - height_error) * pixelWidth; double sx = ((ix - xc) - width_error) * pixelWidth; @@ -619,14 +711,14 @@ void Image::draw_cube(const double *x, const double *surfaceColor, double diamet ------------------------------------------------------------------------- */ void Image::draw_cylinder(const double *x, const double *y, - const double *surfaceColor, double diameter, int sflag) + const double *surfaceColor, double diameter, int sflag, double opacity) { double mid[3],xaxis[3],yaxis[3],zaxis[3]; double camLDir[3], camLRight[3], camLUp[3]; double zmin, zmax; - if (sflag % 2) draw_sphere(x,surfaceColor,diameter); - if (sflag / 2) draw_sphere(y,surfaceColor,diameter); + if (sflag % 2) draw_sphere(x,surfaceColor,diameter,opacity); + if (sflag / 2) draw_sphere(y,surfaceColor,diameter,opacity); double radius = 0.5*diameter; double radsq = radius*radius; @@ -700,6 +792,7 @@ void Image::draw_cylinder(const double *x, const double *y, for (int iy = yc - pixelHalfHeight; iy <= yc + pixelHalfHeight; iy ++) { for (int ix = xc - pixelHalfWidth; ix <= xc + pixelHalfWidth; ix ++) { if (iy < 0 || iy >= height || ix < 0 || ix >= width) continue; + if (((opacity < 1.0) && (transthresh[ix % TRANK][iy % TRANK] > opacity)) || (opacity <= 0.0)) continue; double surface[3], normal[3]; double sy = ((iy - yc) - height_error) * pixelWidth; @@ -744,10 +837,11 @@ void Image::draw_cylinder(const double *x, const double *y, } /* ---------------------------------------------------------------------- - draw triangle with 3 corner points x,y,z and surfaceColor + draw triangle with 3 corner points x,y,z, surfaceColor ------------------------------------------------------------------------- */ -void Image::draw_triangle(const double *x, const double *y, const double *z, const double *surfaceColor) +void Image::draw_triangle(const double *x, const double *y, const double *z, const double *surfaceColor, + const double opacity) { double d1[3], d1len, d2[3], d2len, normal[3], invndotd; double xlocal[3], ylocal[3], zlocal[3]; @@ -825,6 +919,7 @@ void Image::draw_triangle(const double *x, const double *y, const double *z, con for (int iy = yc - pixelDown; iy <= yc + pixelUp; iy ++) { for (int ix = xc - pixelLeft; ix <= xc + pixelRight; ix ++) { if (iy < 0 || iy >= height || ix < 0 || ix >= width) continue; + if (((opacity < 1.0) && (transthresh[ix % TRANK][iy % TRANK] > opacity)) || (opacity <= 0.0)) continue; double sy = ((iy - yc) - height_error) * pixelWidth; double sx = ((ix - xc) - width_error) * pixelWidth; @@ -1557,6 +1652,7 @@ double *Image::color2rgb(const char *color, int index) } if (color) { + if (strcmp(color,"none") == 0) return nullptr; for (int i = 0; i < ncolors; i++) if (strcmp(color,username[i]) == 0) return userrgb[i]; for (int i = 0; i < NCOLORS; i++) diff --git a/src/image.h b/src/image.h index 71d671c0414..eecfe29564b 100644 --- a/src/image.h +++ b/src/image.h @@ -50,12 +50,14 @@ class Image : protected Pointers { void write_PPM(FILE *); void view_params(double, double, double, double, double, double); - void draw_sphere(const double *, const double *, double); - void draw_cube(const double *, const double *, double); - void draw_cylinder(const double *, const double *, const double *, double, int); - void draw_triangle(const double *, const double *, const double *, const double *); - void draw_box(double (*)[3], double); - void draw_axes(double (*)[3], double); + void draw_sphere(const double *, const double *, double, double opacity = 1.0); + void draw_cube(const double *, const double *, double, double opacity = 1.0); + void draw_cylinder(const double *, const double *, const double *, double, int, + double opacity = 1.0); + void draw_triangle(const double *, const double *, const double *, const double *, + double opacity = 1.0); + void draw_box(double (*)[3], double, double opacity = 1.0); + void draw_axes(double (*)[3], double, double opacity = 1.0); int map_dynamic(int); int map_reset(int, int, char **);