Dynamic Geometry

Overview

Before we get into how the algorithm works, here is how it is used:

Mesh "snapshots" are referenced in a text file similar to how variable rate constant files are handled. The first column contains MDL filenames and the second column contains the time at which the geometry should be used. As an example, imagine we have a file called geometry_snapshots.txt that contained the following:

my_geometry0.mdl 0e-6
my_geometry1.mdl 1e-6
my_geometry2.mdl 2e-6
my_geometry3.mdl 3e-6
...

Each MDL listed would contain mesh object definitions (usually one wouldn't add new meshes or remove existing ones, but this is now allowed) where the vertices and elements might be different. You would also need to re-instantiate the existing mesh objects. For example, my_geometry0.mdl might contain the following:

Cube POLYGON_LIST {
  VERTEX_LIST {
    [  0.1,  0.1, -0.1 ]
    ...
  }
  ELEMENT_CONNECTIONS {
    [ 0, 1, 2 ]
    ...
  }
}
INSTANTIATE Scene OBJECT
{
  Cube OBJECT Cube{}
}

In this example, each subsequent file referenced in geometry_snapshots.txt (e.g. my_geometry1.mdl, my_geometry2.mdl, etc) might have a definition and instantiation of the Cube object. As mentioned before, the actual vertices and elements can be different.

Finally, the file containing this list of MDLs and times (geometry_snapshots.txt in this example) would be referenced somewhere in your main MDL using a new keyword called DYNAMIC_GEOMETRY like this:

DYNAMIC_GEOMETRY = "geometry_snapshots.txt"

Note: For lack of a better term, we will refer to the file assigned to DYNAMIC_GEOMETRY file as the "dynamic geometry file". Do not confuse this with the MDLs that have the actual geometry/meshes and instantiations.

The usage is probably best explained in the wiki:

https://github.com/mcellteam/mcell/wiki/Using-Dynamic-Geometries

Algorithm

Now, we can move onto how the algorithm actually works:

Initialization

Along with all the other data structures, memory for dynamic geometry events is allocated in init_data_structures (init.c).

The scheduler is created in init_dynamic_geometry (init.c).

Parsing

During parse time, we check what text file is assigned to DYNAMIC_GEOMETRY (mcell_add_dynamic_geometry_file in mcell_dyngeom.c). We then process the dynamic geometry file, adding an event to the dynamic_geometry_head (add_dynamic_geometry_events in dyngeom.c) for each line in the DG file. Also, we do a preliminary parse of all the geometry files listed in the DG file (parse_dg in dyngeom_yacc.c). We need to do this to make sure that objects we are counting in/on will exist at some point during the simulation. NOTE: This is different than the regular geometry parsing done in mcell_parse_mdl and uses a different parser. The preliminary parsing is somewhat involved and so is covered more thoroughly in a later section of this document. After the preliminary parsing is done, everything in dynamic_geometry_head is added to the scheduler (end of schedule_dynamic_geometry in init.c).

Set Up Some Data Structures

Create a data structure so we can quickly check if a molecule species can move in or out of any given surface region. (init_species_mesh_transp in dyngeom.c)

Run Simulation

Once the simulation starts, check for a scheduled geometry change every timestep (process_geometry_changes in mcell_run.c). If a geometry change is scheduled, we do the following:

Save all the molecules available in the scheduler (save_all_molecules)
- For volume molecules, keep track of enclosing meshes in the order of nesting. (save_volume_molecule)
- For surface molecules, keep track of the mesh it is on as well as all regions. (save_surface_molecule)
- Save common properties like next uni reaction, scheduling time, birthday, etc. (save_common_molecule_properties).
- Side note: Finding what mesh a molecule is inside of is explained in point 2 below. (find_enclosing_meshes)
Save a list of names of all the meshes and regions in fully qualified form prior to trashing them (create_mesh_instantiation_sb). We'll get back to this in a minute.
Now the fun part; destroy geometry, subvolumes, memory helpers, and pretty much the entire simulation (destroy_everything in mcell_redo_geom). We also zero out counts (reset_current_counts), because they'll get repopulated when we count from scratch later. Well, we leave things like molecule and reaction definitions alone, since they're not hurting anyone.
Parse the new geometry and instantiation information. (mcell_parse_mdl in mcell_redo_geom)
Re-initialize the things we just destroyed (init_partitions, init_vertices_walls, init_regions).
Save a list of the new mesh names (create_mesh_instantiation_sb) and compare them with the old names to see if anything has been added or removed (sym_diff_string_buffers). We will ignore newly added or removed meshes and regions for the purpose of placing molecules this time. New meshes will be tracked during subsequent dynamic geometry events.
Place all the molecules that we previously saved, moving molecules if necessary to keep them in/on the appropriate compartment (place_all_molecules).

Placement of Volume Molecules

To expand on the last point, here's how placement works for volume molecules (insert_volume_molecule_encl_mesh):

Find what subvolume the molecule should be in now based on its saved location. (find_subvolume in vol_util.c)
Find the current enclosing meshes. (find_enclosing_meshes in dyngeom.c) A. Pick a random vector and scale it so that it's big enough to cross the entire world. B. Ray trace from the current molecule location along the direction of the random vector. C. Start working through all our collisions a. Keep track of how many times we cross each closed mesh. If odd number, then we are inside it. If even, then we are outside it. (http://en.wikipedia.org/wiki/Point_in_polygon) b. If we reach the edge of the world (outermost subvolume), return each enclosing mesh name (i.e. the first) we found, assuming there is one. If there isn't one, then we are outside all objects and just return NULL. c. If we ever hit an interior subvolume, then repeat 2.B and 2.C until we hit the edge of the world.
Compare where the molecule was (prior to the geometry change) with where it is now relative to the meshes. This also takes into account meshes nested inside other meshes. Ultimately, we want to see if movement is possible from the starting position to the ending position. (compare_molecule_nesting)

If the compartment nesting never changed and there were no molecule transparencies, then we could ignore the nesting and simply track the closest enclosing mesh. In fact, this is what was done in an earlier version of the dynamic geometries code. But that is not the case. In fact, meshes can be selectively transparent to a molecule (in-to-out, out-to-in, both, or neither), the transparent properties can change from one dynamic geometry event to the next, and entire meshes/regions can disappear.

Now, a little discussion on notation for the sake of being clear and concise. When we say something like A->B->C->null, this means that mesh A is inside mesh B which is inside mesh C. Lastly, C is an outermost mesh.

First, find "overlap" if any exists. (check_overlapping_meshes) For example:

 mesh_names_old: A->B->C->D->null
 mesh_names_new:       C->D->null

"A" was the closest enclosing mesh and then "C" became the closest enclosing mesh. That means the molecule moved from "A" to "C".

The order could be reversed like this:

 mesh_names_old:       C->D->null
 mesh_names_new: A->B->C->D->null

This means the molecule moved from "C" to "A".

We could also have a case with no overlap (aside from null) like this:

 mesh_names_old: C->D->null
 mesh_names_new: E->F->null

This means the molecule moved from "C" to "E". (check_nonoverlapping_meshes)

See if we can legally be where we are at (check_outin_or_inout)
If the nesting is the same, just place the molecule where it is.
If the nesting is different, we need to move the molecule (place_mol_relative_to_mesh). A. Go through all the walls in the subvolume that is shared with the molecule. B. Find the distance to the closest wall. C. Check neighboring subvolumes to see if there's an even closer wall. D. Find closest wall and put the molecule either slightly in front of or behind it. By default, this is done to the nearest point (DYNAMIC_GEOMETRY_MOLECULE_PLACEMENT = NEAREST_POINT), but it's also possible to just set it to a random point on the nearest triangle (DYNAMIC_GEOMETRY_MOLECULE_PLACEMENT = NEAREST_TRIANGLE).
Update subvolume, enable counting, schedule it, and do other book keeping.

Placement of Surface Molecules

Now, here's how placement works for surface molecules (insert_surface_molecule and place_surface_molecule in vol_util.c):

Find what subvolume the molecule is now in based on its saved location.
Go through every wall in the subvolume. (find_subvolume) A. Make sure the object name and region names match what we saved previously aside for newly-added/removed regions and regions that we are transparent to. (verify_wall_regions_match in grid_util.c) B. If they match, find the closest point on the wall to our molecule's location (closest_interior_point in wall_util.c). C. Of all the walls, keep the best (i.e. the closest) location.
Make sure there's not a better/closer location in adjacent subvolumes.
Then make sure there's actually a free slot on the wall's surface grid.
If not, check neighboring walls.
If everything is okay, put it on the grid.
Make it countable for reaction output.
Add it to the scheduler.

Preliminary Parsing Geometry

Before parsing, we call create_dg_parse which creates a global (BOO!) pointer to a struct called dg_parse, which stores some useful information about the object(s) we are currently parsing, symbol tables, etc. Next, we call init_top_level_objs, which create the aforementioned symbol tables. These are used to store the objects and regions. We also create two top level objects (WORLD_OBJ and WORLD_INSTANCE), add them to object symbol table, and add them to the top of the root_object and root_instance trees.

The actual parsing itself begins in parse_dg where we set the root object and root instance (setup_root_obj_inst). NOTE: This is unnecessary the first time through, since it's also done in init_top_level_objs. Then we open the geometry files for parsing by calling yyparse.

The first thing we will likely encounter when parsing a geometry file is a polygon_list object (dg_new_polygon_list in dyngeom_parse_extras.c). We prepare the object_name_list. This is necessary in case the object is part of a meta object (this needs tested more thoroughly), such that the fully qualified name is something like Vesicles.Left.Ves1. Then we call dg_start_object_simple (dyngeom_parse_extras.c), which adds the name to the object_name_list to form the fully qualified name and makes the new object (make_new_object). This is where the object name gets added to the DG object symbol table. Next we create an "ALL" region, which every polygon object has (dg_create_region and dg_make_new_region). We make a fully qualified region name of the form "object_name,region_name" and store the name in the DG region symbol table. We may process other regions the user created at this time. When all of this is finished, we call dg_finish_object which basically just pops the name off the object_name_list.

The next major thing that the user is likely to do is instantiate the objects. As every object is created we first call dg_start_object. This is very similar to dg_start_object_simple. We add the object name to the object_name_list and store it in the DG object symbol table.

Eventually, we will copy the contents of one object into another one using dg_deep_copy_object. For normal cases, this mainly entails copying the regions of the source/original object to the destination/new object (dg_copy_object_regions).

dynamic_geometry.md 12 KB 文件历史 原始文件