Example 11 with MembraneSubDomain
use of cbit.vcell.math.MembraneSubDomain in project vcell by virtualcell.
the class SmoldynFileWriter method writeDrifts.
private void writeDrifts() throws ExpressionBindingException, ExpressionException, MathException {
// writer diffusion properties
printWriter.println("# drift properties");
Enumeration<SubDomain> subDomainEnumeration = mathDesc.getSubDomains();
while (subDomainEnumeration.hasMoreElements()) {
SubDomain subDomain = subDomainEnumeration.nextElement();
List<ParticleProperties> particlePropertiesList = subDomain.getParticleProperties();
for (ParticleProperties pp : particlePropertiesList) {
String variableName = null;
if (subDomain instanceof MembraneSubDomain) {
variableName = getVariableName(pp.getVariable(), subDomain);
} else {
variableName = getVariableName(pp.getVariable(), null);
}
try {
double driftX = 0.0;
if (pp.getDriftX() != null) {
driftX = subsituteFlattenToConstant(pp.getDriftX());
}
double driftY = 0.0;
if (pp.getDriftY() != null) {
driftY = subsituteFlattenToConstant(pp.getDriftY());
}
double driftZ = 0.0;
if (pp.getDriftZ() != null) {
driftZ = subsituteFlattenToConstant(pp.getDriftZ());
}
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.drift + " " + variableName + " " + driftX + " " + driftY + " " + driftZ);
} catch (NotAConstantException ex) {
throw new ExpressionException("diffusion coefficient for variable " + variableName + " is not a constant. Constants are required for all diffusion coefficients");
}
}
}
printWriter.println();
}
Also used :
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
ParticleProperties(cbit.vcell.math.ParticleProperties)
ExpressionException(cbit.vcell.parser.ExpressionException)
Example 12 with MembraneSubDomain
use of cbit.vcell.math.MembraneSubDomain in project vcell by virtualcell.
the class SmoldynFileWriter method writeReactions.
private void writeReactions() throws ExpressionException, MathException {
printWriter.println("# reactions");
Enumeration<SubDomain> subdomains = mathDesc.getSubDomains();
while (subdomains.hasMoreElements()) {
SubDomain subdomain = subdomains.nextElement();
for (ParticleJumpProcess pjp : subdomain.getParticleJumpProcesses()) {
ArrayList<Variable> reactants = new ArrayList<Variable>();
ArrayList<Variable> products = new ArrayList<Variable>();
for (Action a : pjp.getActions().toArray(new Action[pjp.getActions().size()])) {
if (a.getOperation().equals(Action.ACTION_CREATE)) {
products.add(a.getVar());
} else if (a.getOperation().equals(Action.ACTION_DESTROY)) {
reactants.add(a.getVar());
}
}
Expression rateDefinition = null;
JumpProcessRateDefinition jprd = pjp.getParticleRateDefinition();
if (jprd instanceof MacroscopicRateConstant) {
rateDefinition = subsituteFlatten(((MacroscopicRateConstant) jprd).getExpression());
} else if (jprd instanceof InteractionRadius) {
rateDefinition = subsituteFlatten(((InteractionRadius) jprd).getExpression());
} else {
new RuntimeException("The jump process rate definition is not supported");
}
if (rateDefinition.isZero()) {
continue;
}
if (mathDesc.isSpatialHybrid()) {
String[] symbols = rateDefinition.getSymbols();
if (symbols != null) {
if (subdomain instanceof MembraneSubDomain) {
rateDefinition = new Expression(FiniteVolumeFileWriter.replaceVolumeVariable(getSimulationTask(), (MembraneSubDomain) subdomain, rateDefinition));
}
}
} else {
try {
rateDefinition.evaluateConstant();
} catch (ExpressionException ex) {
throw new ExpressionException("reaction rate for jump process " + pjp.getName() + " is not a constant. Constants are required for all reaction rates.");
}
}
// Smoldyn takes maximum 2nd order reaction.
if (reactants.size() > 2) {
throw new MathException("VCell spatial stochastic models support up to 2nd order reactions. \n" + "The reaction:" + pjp.getName() + " has more than 2 reactants.");
}
if (products.size() > 2) {
throw new MathException("VCell spatial stochastic models support up to 2nd order reactions. \n" + "The reaction:" + pjp.getName() + " has more than 2 products.");
}
String rateDefinitionStr = simulation.getMathDescription().isSpatialHybrid() ? rateDefinition.infix() + ";" : rateDefinition.evaluateConstant() + "";
if (subdomain instanceof CompartmentSubDomain) {
// 0th order reaction, product limited to one and we'll let the reaction know where it happens
if (reactants.size() == 0 && products.size() == 1) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_cmpt + " " + subdomain.getName() + " " + pjp.getName() + " ");
} else {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction + " " + /* + subdomain.getName() + " "*/
pjp.getName() + " ");
}
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else if (subdomain instanceof MembraneSubDomain) {
// 0th order reaction, product limited to one and it can be on mem or in vol
if (reactants.size() == 0 && products.size() == 1) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else // consuming of a species to nothing, limited to one reactant
if (reactants.size() == 1 && products.size() == 0) {
if (// consuming a mem species in mem reaction
getMembraneVariableCount(reactants) == 1) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else // it equals to adsorption, species A from vol adsorbed to mem as again species A, and then we kill the speceis A on mem.
if (getVolumeVariableCount(reactants) == 1) {
writeRateTransitionCommand(reactants, products, subdomain, rateDefinitionStr);
String speciesName = reactants.get(0).getName();
String killMolCmd = "cmd " + SmoldynVCellMapper.SmoldynKeyword.E + " " + SmoldynVCellMapper.SmoldynKeyword.killmol + " " + speciesName + "(" + SmoldynKeyword.up + ")";
killMolCommands.add(killMolCmd);
}
} else // Use rate command for any membrane reactions with 1 reactant and 1 product
if ((reactants.size() == 1) && (products.size() == 1)) {
// Membrane reaction (1 react to 1 product).
if (getMembraneVariableCount(products) == 1 && getMembraneVariableCount(reactants) == 1) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else // Other single molecular reactions
{
writeRateTransitionCommand(reactants, products, subdomain, rateDefinitionStr);
}
} else // membrane reactions which are not one to one, or 0th order, or consuming species
{
if (// membrane reaction has one membrane bound reactant
(getMembraneVariableCount(reactants) == 1)) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else if (// bimolecular membrane reaction
getMembraneVariableCount(reactants) == 2) {
if (jprd instanceof InteractionRadius) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionByInteractionRadius(reactants, products, subdomain, rateDefinitionStr, pjp.getName());
} else {
// throw new MathException("Error with reaction: " + pjp.getName() + ".\nVCell Spatial stochastic modeling requires macroscopic or microscopic kinetics for bimolecular membrane reactions.");
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
}
} else if (getMembraneVariableCount(reactants) == 0) {
throw new MathException("Error with reaction: " + pjp.getName() + ".\nIn VCell spatial stochastic modeling, the membrane reaction requires at least one membrane bound reactant.");
}
}
}
}
}
printWriter.println();
}
Also used :
Action(cbit.vcell.math.Action)
JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
ReservedVariable(cbit.vcell.math.ReservedVariable)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
Variable(cbit.vcell.math.Variable)
InteractionRadius(cbit.vcell.math.InteractionRadius)
ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess)
ArrayList(java.util.ArrayList)
ExpressionException(cbit.vcell.parser.ExpressionException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
Expression(cbit.vcell.parser.Expression)
MathException(cbit.vcell.math.MathException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
Example 13 with MembraneSubDomain
use of cbit.vcell.math.MembraneSubDomain in project vcell by virtualcell.
the class SmoldynSurfaceTessellator method writeSurfaces.
protected void writeSurfaces() throws SolverException, ImageException, PropertyVetoException, GeometryException, ExpressionException {
GeometrySurfaceDescription geometrySurfaceDescription = resampledGeometry.getGeometrySurfaceDescription();
SurfaceClass[] surfaceClasses = geometrySurfaceDescription.getSurfaceClasses();
GeometrySpec geometrySpec = resampledGeometry.getGeometrySpec();
SubVolume[] surfaceGeometrySubVolumes = geometrySpec.getSubVolumes();
GeometricRegion[] AllGeometricRegions = resampledGeometry.getGeometrySurfaceDescription().getGeometricRegions();
ArrayList<SurfaceGeometricRegion> surfaceRegionList = new ArrayList<SurfaceGeometricRegion>();
ArrayList<VolumeGeometricRegion> volumeRegionList = new ArrayList<VolumeGeometricRegion>();
for (GeometricRegion geometricRegion : AllGeometricRegions) {
if (geometricRegion instanceof SurfaceGeometricRegion) {
surfaceRegionList.add((SurfaceGeometricRegion) geometricRegion);
} else if (geometricRegion instanceof VolumeGeometricRegion) {
volumeRegionList.add((VolumeGeometricRegion) geometricRegion);
} else {
throw new SolverException("unsupported geometric region type " + geometricRegion.getClass());
}
}
printWriter.println("# geometry");
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.dim + " " + dimension);
if (bHasNoSurface) {
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_compartment + " " + surfaceGeometrySubVolumes.length);
} else {
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_compartment + " " + (surfaceGeometrySubVolumes.length + 1));
// plus the surface which are bounding walls
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_surface + " " + (surfaceClasses.length + dimension));
}
printWriter.println();
// write boundaries and wall surfaces
writeWallSurfaces();
// for 3D ... smoldyn normal convension is triangle right-hand-rule normal points to the outside compartment subdomain.
if (!bHasNoSurface) {
membraneSubdomainTriangleMap = new HashMap<MembraneSubDomain, ArrayList<TrianglePanel>>();
// write surfaces
printWriter.println("# surfaces");
int triangleGlobalCount = 0;
int membraneIndex = -1;
SurfaceCollection surfaceCollection = geometrySurfaceDescription.getSurfaceCollection();
// pre-allocate collections used repeatedly in following loops; clear before reusing
HashMap<Node, Set<String>> nodeTriMap = new HashMap<>();
ArrayList<TrianglePanel> triList = new ArrayList<TrianglePanel>();
// use a sorted set to ensure neighbors written out is same order for reproducibility
SortedSet<String> neighborsForCurrentNode = new TreeSet<String>();
for (int sci = 0; sci < surfaceClasses.length; sci++) {
nodeTriMap.clear();
triList.clear();
int triLocalCount = 0;
SurfaceClass surfaceClass = surfaceClasses[sci];
GeometricRegion[] geometricRegions = geometrySurfaceDescription.getGeometricRegions(surfaceClass);
for (GeometricRegion gr : geometricRegions) {
SurfaceGeometricRegion sgr = (SurfaceGeometricRegion) gr;
VolumeGeometricRegion volRegion0 = (VolumeGeometricRegion) sgr.getAdjacentGeometricRegions()[0];
VolumeGeometricRegion volRegion1 = (VolumeGeometricRegion) sgr.getAdjacentGeometricRegions()[1];
SubVolume subVolume0 = volRegion0.getSubVolume();
SubVolume subVolume1 = volRegion1.getSubVolume();
CompartmentSubDomain compart0 = mathDesc.getCompartmentSubDomain(subVolume0.getName());
CompartmentSubDomain compart1 = mathDesc.getCompartmentSubDomain(subVolume1.getName());
MembraneSubDomain membraneSubDomain = mathDesc.getMembraneSubDomain(compart0, compart1);
if (membraneSubDomain == null) {
throw new SolverException(VCellErrorMessages.getSmoldynUnexpectedSurface(compart0, compart1));
}
int exteriorRegionID = volRegion0.getRegionID();
int interiorRegionID = volRegion1.getRegionID();
if (membraneSubDomain.getInsideCompartment() == compart0) {
exteriorRegionID = volRegion1.getRegionID();
interiorRegionID = volRegion0.getRegionID();
}
for (int j = 0; j < surfaceCollection.getSurfaceCount(); j++) {
Surface surface = surfaceCollection.getSurfaces(j);
if ((surface.getInteriorRegionIndex() == exteriorRegionID && surface.getExteriorRegionIndex() == interiorRegionID) || (surface.getInteriorRegionIndex() == interiorRegionID && surface.getExteriorRegionIndex() == exteriorRegionID)) {
// Polygon polygon = surface.getPolygons(k);
for (Polygon polygon : surface) {
if (polygonMembaneElementMap != null) {
membraneIndex = polygonMembaneElementMap.get(polygon).getMembraneIndex();
}
Node[] nodes = polygon.getNodes();
if (dimension == 2) {
// ignore z
Vect3d unitNormal = new Vect3d();
polygon.getUnitNormal(unitNormal);
unitNormal.set(unitNormal.getX(), unitNormal.getY(), 0);
int point0 = 0;
Vect3d v0 = new Vect3d(nodes[point0].getX(), nodes[point0].getY(), 0);
int point1 = 1;
Vect3d v1 = null;
for (point1 = 1; point1 < nodes.length; point1++) {
if (v0.getX() != nodes[point1].getX() || v0.getY() != nodes[point1].getY()) {
v1 = new Vect3d(nodes[point1].getX(), nodes[point1].getY(), 0);
break;
}
}
if (v1 == null) {
throw new RuntimeException("failed to generate surface");
}
Vect3d v01 = Vect3d.sub(v1, v0);
Vect3d unit01n = v01.cross(unitNormal);
unit01n.unit();
if (Math.abs(unit01n.getZ() - 1.0) < 1e-6) {
// v0 to v1 opposes vcell surface normal. it's already flipped.
Triangle triangle;
if (surface.getInteriorRegionIndex() == interiorRegionID) {
// we have to flipped it back
triangle = new Triangle(nodes[point1], nodes[point0], null);
} else {
triangle = new Triangle(nodes[point0], nodes[point1], null);
}
triList.add(new TrianglePanel(triLocalCount++, triangleGlobalCount++, membraneIndex, triangle));
} else if (Math.abs(unit01n.getZ() + 1.0) < 1e-6) {
// v0 to v1 is in direction of vcell surface normal.
Triangle triangle;
if (surface.getInteriorRegionIndex() == interiorRegionID) {
triangle = new Triangle(nodes[point0], nodes[point1], null);
} else {
triangle = new Triangle(nodes[point1], nodes[point0], null);
}
triList.add(new TrianglePanel(triLocalCount++, triangleGlobalCount++, membraneIndex, triangle));
} else {
throw new RuntimeException("failed to generate surface");
}
} else if (dimension == 3) {
Triangle triangle1;
Triangle triangle2;
if (surface.getInteriorRegionIndex() == interiorRegionID) {
// interior
triangle1 = new Triangle(nodes[0], nodes[1], nodes[2]);
triangle2 = new Triangle(nodes[0], nodes[2], nodes[3]);
} else {
triangle1 = new Triangle(nodes[2], nodes[1], nodes[0]);
triangle2 = new Triangle(nodes[3], nodes[2], nodes[0]);
}
triList.add(new TrianglePanel(triLocalCount++, triangleGlobalCount++, membraneIndex, triangle1));
triList.add(new TrianglePanel(triLocalCount++, triangleGlobalCount++, membraneIndex, triangle2));
}
}
}
}
}
// add triangles to node hash
for (TrianglePanel triPanel : triList) {
for (Node node : triPanel.triangle.getNodes()) {
if (node == null) {
continue;
}
Set<String> triNameSet = nodeTriMap.get(node);
if (triNameSet == null) {
triNameSet = new HashSet<String>();
nodeTriMap.put(node, triNameSet);
}
triNameSet.add(triPanel.name);
}
}
SubVolume[] adjacentSubvolums = surfaceClass.getAdjacentSubvolumes().toArray(new SubVolume[0]);
CompartmentSubDomain csd0 = simulation.getMathDescription().getCompartmentSubDomain(adjacentSubvolums[0].getName());
CompartmentSubDomain csd1 = simulation.getMathDescription().getCompartmentSubDomain(adjacentSubvolums[1].getName());
MembraneSubDomain membraneSubDomain = simulation.getMathDescription().getMembraneSubDomain(csd0, csd1);
membraneSubdomainTriangleMap.put(membraneSubDomain, triList);
final boolean initialMoleculesOnMembrane = (closestTriangles != null);
if (initialMoleculesOnMembrane) {
findClosestTriangles(membraneSubDomain, triList);
}
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.start_surface + " " + surfaceClass.getName());
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + "(" + SmoldynVCellMapper.SmoldynKeyword.all + ") " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.reflect);
// printWriter.println(SmoldynKeyword.action + " " + SmoldynKeyword.all + "(" + SmoldynKeyword.up + ") " + SmoldynKeyword.both + " " + SmoldynKeyword.reflect);
Color c = colorForSurface(sci);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.color + " " + SmoldynVCellMapper.SmoldynKeyword.both + " " + c.getRed() / 255.0 + " " + c.getGreen() / 255.0 + " " + c.getBlue() / 255.0 + " 0.1");
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.front + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.back + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_panels + " " + SmoldynVCellMapper.SmoldynKeyword.tri + " " + triList.size());
for (TrianglePanel trianglePanel : triList) {
Triangle triangle = trianglePanel.triangle;
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.tri);
switch(dimension) {
case 1:
printWriter.print(" " + triangle.getNodes(0).getX());
break;
case 2:
printWriter.print(" " + triangle.getNodes(0).getX() + " " + triangle.getNodes(0).getY());
printWriter.print(" " + triangle.getNodes(1).getX() + " " + triangle.getNodes(1).getY());
break;
case 3:
for (Node node : triangle.getNodes()) {
printWriter.print(" " + node.getX() + " " + node.getY() + " " + node.getZ());
}
break;
}
printWriter.println(" " + trianglePanel.name);
}
for (TrianglePanel triPanel : triList) {
neighborsForCurrentNode.clear();
for (Node node : triPanel.triangle.getNodes()) {
if (node == null) {
continue;
}
neighborsForCurrentNode.addAll(nodeTriMap.get(node));
}
neighborsForCurrentNode.remove(triPanel.name);
// printWriter.print(SmoldynKeyword.neighbors + " " +triPanel.name);
// to allow smoldyn read line length as 256, chop the neighbors to multiple lines
int maxNeighborCount = 4;
//
int count = 0;
for (String neigh : neighborsForCurrentNode) {
if (count % maxNeighborCount == 0) {
printWriter.println();
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.neighbors + " " + triPanel.name);
}
printWriter.print(" " + neigh);
count++;
}
}
printWriter.println();
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.end_surface);
printWriter.println();
}
}
}
Also used :
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
SortedSet(java.util.SortedSet)
TreeSet(java.util.TreeSet)
HashSet(java.util.HashSet)
Set(java.util.Set)
GeometrySurfaceDescription(cbit.vcell.geometry.surface.GeometrySurfaceDescription)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
HashMap(java.util.HashMap)
Node(cbit.vcell.geometry.surface.Node)
ArrayList(java.util.ArrayList)
Triangle(cbit.vcell.geometry.surface.Triangle)
Surface(cbit.vcell.geometry.surface.Surface)
GeometrySpec(cbit.vcell.geometry.GeometrySpec)
SubVolume(cbit.vcell.geometry.SubVolume)
TreeSet(java.util.TreeSet)
Polygon(cbit.vcell.geometry.surface.Polygon)
SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion)
SurfaceCollection(cbit.vcell.geometry.surface.SurfaceCollection)
Color(java.awt.Color)
VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion)
VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion)
SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion)
GeometricRegion(cbit.vcell.geometry.surface.GeometricRegion)
Vect3d(cbit.vcell.render.Vect3d)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SolverException(cbit.vcell.solver.SolverException)
Example 14 with MembraneSubDomain
use of cbit.vcell.math.MembraneSubDomain in project vcell by virtualcell.
the class MathTestingUtilities method comparePDEResultsWithExact.
/**
* Insert the method's description here.
* Creation date: (8/20/2003 12:58:10 PM)
*/
public static SimulationComparisonSummary comparePDEResultsWithExact(SimulationSymbolTable simSymbolTable, PDEDataManager dataManager, String type, double absErrorThreshold, double relErrorThreshold) throws DataAccessException, ExpressionException {
java.util.Hashtable<String, DataErrorSummary> tempVarHash = new java.util.Hashtable<String, DataErrorSummary>();
double[] timeArray = dataManager.getDataSetTimes();
Variable[] vars = simSymbolTable.getVariables();
CartesianMesh mesh = dataManager.getMesh();
MathDescription mathDesc = simSymbolTable.getSimulation().getMathDescription();
// Get volumeSubdomains from mathDesc/mesh and store in lookupTable
int numVol = mesh.getSizeX() * mesh.getSizeY() * mesh.getSizeZ();
CompartmentSubDomain[] volSubDomainLookup = new CompartmentSubDomain[numVol];
for (int i = 0; i < numVol; i++) {
int subVolumeIndex = mesh.getSubVolumeFromVolumeIndex(i);
SubVolume subVolume = mathDesc.getGeometry().getGeometrySpec().getSubVolume(subVolumeIndex);
CompartmentSubDomain compSubDomain = mathDesc.getCompartmentSubDomain(subVolume.getName());
volSubDomainLookup[i] = compSubDomain;
}
// Get membraneSubdomains from mathDesc/mesh and store in lookupTable
int numMem = mesh.getMembraneElements().length;
MembraneSubDomain[] memSubDomainLookup = new MembraneSubDomain[numMem];
for (int i = 0; i < numMem; i++) {
int insideVolIndex = mesh.getMembraneElements()[i].getInsideVolumeIndex();
int outsideVolIndex = mesh.getMembraneElements()[i].getOutsideVolumeIndex();
MembraneSubDomain memSubDomain = mathDesc.getMembraneSubDomain(volSubDomainLookup[insideVolIndex], volSubDomainLookup[outsideVolIndex]);
memSubDomainLookup[i] = memSubDomain;
}
double[] valueArray = new double[4];
SimpleSymbolTable symbolTable = new SimpleSymbolTable(new String[] { "t", "x", "y", "z" });
int tIndex = symbolTable.getEntry("t").getIndex();
int xIndex = symbolTable.getEntry("x").getIndex();
int yIndex = symbolTable.getEntry("y").getIndex();
int zIndex = symbolTable.getEntry("z").getIndex();
SimulationComparisonSummary simComparisonSummary = new SimulationComparisonSummary();
String hashKey = new String("");
long dataLength = 0;
// for each var, do the following :
for (int i = 0; i < vars.length; i++) {
if (vars[i] instanceof VolVariable || vars[i] instanceof MemVariable || vars[i] instanceof FilamentVariable || vars[i] instanceof VolumeRegionVariable || vars[i] instanceof MembraneRegionVariable || vars[i] instanceof FilamentRegionVariable) {
// for each time in timeArray,
for (int j = 0; j < timeArray.length; j++) {
if (type.equals(TestCaseNew.EXACT_STEADY)) {
if (j != (timeArray.length - 1)) {
continue;
}
}
// get data block from varName, data from datablock
SimDataBlock simDataBlock = dataManager.getSimDataBlock(vars[i].getName(), timeArray[j]);
double[] data = simDataBlock.getData();
dataLength = data.length;
SubDomain subDomain = null;
Coordinate subDomainCoord = null;
// for each point in data block ...
for (int k = 0; k < dataLength; k++) {
// Get subdomain from mesh (from the lookupTable), get coordinates (x,y,z) from mesh, evaluate EXACT SOLN at that coord
if (vars[i] instanceof VolVariable) {
subDomain = volSubDomainLookup[k];
subDomainCoord = mesh.getCoordinateFromVolumeIndex(k);
} else if (vars[i] instanceof MemVariable) {
subDomain = memSubDomainLookup[k];
subDomainCoord = mesh.getCoordinateFromMembraneIndex(k);
} else {
throw new RuntimeException("Var " + vars[i].getName() + " not supported yet!");
}
hashKey = vars[i].getName() + ":" + subDomain.getName();
DataErrorSummary tempVar = (DataErrorSummary) tempVarHash.get(hashKey);
if (tempVar == null) {
Expression exp = new Expression(subDomain.getEquation(vars[i]).getExactSolution());
exp.bindExpression(simSymbolTable);
exp = MathUtilities.substituteFunctions(exp, simSymbolTable);
exp = exp.flatten();
exp.bindExpression(symbolTable);
tempVar = new DataErrorSummary(exp);
tempVarHash.put(hashKey, tempVar);
}
// time
valueArray[tIndex] = timeArray[j];
// x
valueArray[xIndex] = subDomainCoord.getX();
// y
valueArray[yIndex] = subDomainCoord.getY();
// z
valueArray[zIndex] = subDomainCoord.getZ();
// EXACT soln at coord subDomainCoord
double value = tempVar.getExactExp().evaluateVector(valueArray);
tempVar.addDataValues(value, data[k], timeArray[j], k, absErrorThreshold, relErrorThreshold);
}
// end for (k)
}
// end for (j)
}
// end - if (var)
}
// end for (i)
Enumeration<String> enumKeys = tempVarHash.keys();
while (enumKeys.hasMoreElements()) {
String key = enumKeys.nextElement();
DataErrorSummary tempVarSummary = tempVarHash.get(key);
simComparisonSummary.addVariableComparisonSummary(new VariableComparisonSummary(key, tempVarSummary.getMinRef(), tempVarSummary.getMaxRef(), tempVarSummary.getMaxAbsoluteError(), tempVarSummary.getMaxRelativeError(), tempVarSummary.getL2Norm(), tempVarSummary.getTimeAtMaxAbsoluteError(), tempVarSummary.getIndexAtMaxAbsoluteError(), tempVarSummary.getTimeAtMaxRelativeError(), tempVarSummary.getIndexAtMaxRelativeError()));
}
return simComparisonSummary;
}
Also used :
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable)
InsideVariable(cbit.vcell.math.InsideVariable)
SensVariable(cbit.vcell.solver.ode.SensVariable)
FilamentVariable(cbit.vcell.math.FilamentVariable)
VolVariable(cbit.vcell.math.VolVariable)
VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable)
ReservedVariable(cbit.vcell.math.ReservedVariable)
MemVariable(cbit.vcell.math.MemVariable)
OutsideVariable(cbit.vcell.math.OutsideVariable)
FilamentRegionVariable(cbit.vcell.math.FilamentRegionVariable)
Variable(cbit.vcell.math.Variable)
MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable)
MathDescription(cbit.vcell.math.MathDescription)
VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
MemVariable(cbit.vcell.math.MemVariable)
SimDataBlock(cbit.vcell.simdata.SimDataBlock)
SubVolume(cbit.vcell.geometry.SubVolume)
AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume)
VolVariable(cbit.vcell.math.VolVariable)
FilamentRegionVariable(cbit.vcell.math.FilamentRegionVariable)
CartesianMesh(cbit.vcell.solvers.CartesianMesh)
SimpleSymbolTable(cbit.vcell.parser.SimpleSymbolTable)
Coordinate(org.vcell.util.Coordinate)
Expression(cbit.vcell.parser.Expression)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
FilamentVariable(cbit.vcell.math.FilamentVariable)
Example 15 with MembraneSubDomain
use of cbit.vcell.math.MembraneSubDomain in project vcell by virtualcell.
the class MathTestingUtilities method constructExactMath.
/**
* constructExactMath()
*
* take an equation of the form:
*
* d A
* --- = F(A,t)
* d t
*
* and create a new equation with a known exact solution 'A_exact' by adding a forcing function R(t)
*
* d A
* --- = F(A,t) + R(t)
* d t
*
* where:
*
* d A_exact
* R(t) = --------- - F(A_exact,t)
* d t
*
* solving for R(t) is done analytically.
*
* Creation date: (1/21/2003 10:47:54 AM)
* @return cbit.vcell.math.MathDescription
* @param mathDesc cbit.vcell.math.MathDescription
*/
public static MathDescription constructExactMath(MathDescription mathDesc, java.util.Random random, ConstructedSolutionTemplate constructedSolutionTemplate) throws ExpressionException, MathException, MappingException {
if (mathDesc.hasFastSystems()) {
throw new RuntimeException("SolverTest.constructExactMath() suppport for fastSystems not yet implemented.");
}
MathDescription exactMath = null;
try {
exactMath = (MathDescription) BeanUtils.cloneSerializable(mathDesc);
exactMath.setDescription("constructed exact solution from MathDescription (" + mathDesc.getName() + ")");
exactMath.setName("exact from " + mathDesc.getName());
} catch (Throwable e) {
e.printStackTrace(System.out);
throw new RuntimeException("error cloning MathDescription: " + e.getMessage());
}
//
// preload the VariableHash with existing Variables (and Constants,Functions,etc) and then sort all at once.
//
VariableHash varHash = new VariableHash();
Enumeration<Variable> enumVar = exactMath.getVariables();
while (enumVar.hasMoreElements()) {
varHash.addVariable(enumVar.nextElement());
}
java.util.Enumeration<SubDomain> subDomainEnum = exactMath.getSubDomains();
while (subDomainEnum.hasMoreElements()) {
SubDomain subDomain = subDomainEnum.nextElement();
Domain domain = new Domain(subDomain);
java.util.Enumeration<Equation> equationEnum = subDomain.getEquations();
if (subDomain instanceof MembraneSubDomain) {
MembraneSubDomain memSubDomain = (MembraneSubDomain) subDomain;
AnalyticSubVolume insideAnalyticSubVolume = (AnalyticSubVolume) exactMath.getGeometry().getGeometrySpec().getSubVolume(memSubDomain.getInsideCompartment().getName());
Function[] outwardNormalFunctions = getOutwardNormal(insideAnalyticSubVolume.getExpression(), "_" + insideAnalyticSubVolume.getName());
for (int i = 0; i < outwardNormalFunctions.length; i++) {
varHash.addVariable(outwardNormalFunctions[i]);
}
}
while (equationEnum.hasMoreElements()) {
Equation equation = equationEnum.nextElement();
if (equation.getExactSolution() != null) {
throw new RuntimeException("exact solution already exists");
}
Enumeration<Constant> origMathConstants = mathDesc.getConstants();
if (equation instanceof OdeEquation) {
OdeEquation odeEquation = (OdeEquation) equation;
Expression substitutedRateExp = substituteWithExactSolution(odeEquation.getRateExpression(), (CompartmentSubDomain) subDomain, exactMath);
SolutionTemplate solutionTemplate = constructedSolutionTemplate.getSolutionTemplate(equation.getVariable().getName(), subDomain.getName());
String varName = odeEquation.getVariable().getName();
String initName = null;
while (origMathConstants.hasMoreElements()) {
Constant constant = origMathConstants.nextElement();
if (constant.getName().startsWith(varName + "_" + subDomain.getName() + DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONC_UNIT_PREFIX)) {
initName = constant.getName();
}
}
String exactName = varName + "_" + subDomain.getName() + "_exact";
String errorName = varName + "_" + subDomain.getName() + "_error";
String origRateName = "_" + varName + "_" + subDomain.getName() + "_origRate";
String substitutedRateName = "_" + varName + "_" + subDomain.getName() + "_substitutedRate";
String exactTimeDerivativeName = "_" + varName + "_" + subDomain.getName() + "_exact_dt";
Expression exactExp = solutionTemplate.getTemplateExpression();
Expression errorExp = new Expression(exactName + " - " + varName);
Expression origRateExp = new Expression(odeEquation.getRateExpression());
Expression exactTimeDerivativeExp = exactExp.differentiate("t").flatten();
Expression newRate = new Expression(origRateName + " - " + substitutedRateName + " + " + exactTimeDerivativeName);
Constant[] constants = solutionTemplate.getConstants();
for (int i = 0; i < constants.length; i++) {
varHash.addVariable(constants[i]);
}
Expression initExp = new Expression(exactExp);
initExp.substituteInPlace(new Expression("t"), new Expression(0.0));
varHash.addVariable(new Function(initName, initExp.flatten(), domain));
varHash.addVariable(new Function(exactName, exactExp, domain));
varHash.addVariable(new Function(errorName, errorExp, domain));
varHash.addVariable(new Function(exactTimeDerivativeName, exactTimeDerivativeExp, domain));
varHash.addVariable(new Function(origRateName, origRateExp, domain));
varHash.addVariable(new Function(substitutedRateName, substitutedRateExp, domain));
odeEquation.setRateExpression(newRate);
odeEquation.setInitialExpression(new Expression(initName));
odeEquation.setExactSolution(new Expression(exactName));
} else if (equation instanceof PdeEquation) {
PdeEquation pdeEquation = (PdeEquation) equation;
Expression substitutedRateExp = substituteWithExactSolution(pdeEquation.getRateExpression(), (CompartmentSubDomain) subDomain, exactMath);
SolutionTemplate solutionTemplate = constructedSolutionTemplate.getSolutionTemplate(equation.getVariable().getName(), subDomain.getName());
String varName = pdeEquation.getVariable().getName();
String initName = null;
while (origMathConstants.hasMoreElements()) {
Constant constant = origMathConstants.nextElement();
if (constant.getName().startsWith(varName + "_" + subDomain.getName() + DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONC_UNIT_PREFIX)) {
initName = constant.getName();
}
}
String diffusionRateName = "_" + varName + "_" + subDomain.getName() + "_diffusionRate";
String exactName = varName + "_" + subDomain.getName() + "_exact";
String errorName = varName + "_" + subDomain.getName() + "_error";
String origRateName = "_" + varName + "_" + subDomain.getName() + "_origRate";
String substitutedRateName = "_" + varName + "_" + subDomain.getName() + "_substitutedRate";
String exactTimeDerivativeName = "_" + varName + "_" + subDomain.getName() + "_exact_dt";
String exactDxName = "_" + varName + "_" + subDomain.getName() + "_exact_dx";
String exactDyName = "_" + varName + "_" + subDomain.getName() + "_exact_dy";
String exactDzName = "_" + varName + "_" + subDomain.getName() + "_exact_dz";
String exactDx2Name = "_" + varName + "_" + subDomain.getName() + "_exact_dx2";
String exactDy2Name = "_" + varName + "_" + subDomain.getName() + "_exact_dy2";
String exactDz2Name = "_" + varName + "_" + subDomain.getName() + "_exact_dz2";
String exactLaplacianName = "_" + varName + "_" + subDomain.getName() + "_exact_laplacian";
Expression exactExp = solutionTemplate.getTemplateExpression();
Expression errorExp = new Expression(exactName + " - " + varName);
Expression origRateExp = new Expression(pdeEquation.getRateExpression());
Expression initExp = new Expression(exactExp);
initExp.substituteInPlace(new Expression("t"), new Expression(0.0));
initExp = initExp.flatten();
Expression exactTimeDerivativeExp = exactExp.differentiate("t").flatten();
Expression exactDxExp = exactExp.differentiate("x").flatten();
Expression exactDx2Exp = exactDxExp.differentiate("x").flatten();
Expression exactDyExp = exactExp.differentiate("y").flatten();
Expression exactDy2Exp = exactDxExp.differentiate("y").flatten();
Expression exactDzExp = exactExp.differentiate("z").flatten();
Expression exactDz2Exp = exactDxExp.differentiate("z").flatten();
Expression exactLaplacianExp = Expression.add(Expression.add(exactDx2Exp, exactDy2Exp), exactDz2Exp).flatten();
Expression newRate = new Expression(origRateName + " - " + substitutedRateName + " - ((" + diffusionRateName + ")*" + exactLaplacianName + ")" + " + " + exactTimeDerivativeName);
Constant[] constants = solutionTemplate.getConstants();
for (int i = 0; i < constants.length; i++) {
varHash.addVariable(constants[i]);
}
varHash.addVariable(new Function(initName, initExp, domain));
varHash.addVariable(new Function(diffusionRateName, new Expression(pdeEquation.getDiffusionExpression()), domain));
varHash.addVariable(new Function(exactName, exactExp, domain));
varHash.addVariable(new Function(errorName, errorExp, domain));
varHash.addVariable(new Function(exactTimeDerivativeName, exactTimeDerivativeExp, domain));
varHash.addVariable(new Function(origRateName, origRateExp, domain));
varHash.addVariable(new Function(substitutedRateName, substitutedRateExp, domain));
varHash.addVariable(new Function(exactDxName, exactDxExp, domain));
varHash.addVariable(new Function(exactDyName, exactDyExp, domain));
varHash.addVariable(new Function(exactDzName, exactDzExp, domain));
varHash.addVariable(new Function(exactDx2Name, exactDx2Exp, domain));
varHash.addVariable(new Function(exactDy2Name, exactDy2Exp, domain));
varHash.addVariable(new Function(exactDz2Name, exactDz2Exp, domain));
varHash.addVariable(new Function(exactLaplacianName, exactLaplacianExp, domain));
pdeEquation.setRateExpression(newRate);
pdeEquation.setInitialExpression(new Expression(initName));
pdeEquation.setDiffusionExpression(new Expression(diffusionRateName));
pdeEquation.setExactSolution(new Expression(exactName));
CompartmentSubDomain compartmentSubDomain = (CompartmentSubDomain) subDomain;
if (compartmentSubDomain.getBoundaryConditionXm().isDIRICHLET()) {
Expression origExp = pdeEquation.getBoundaryXm();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryXm(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
} else {
pdeEquation.setBoundaryXm(exactExp);
}
} else if (compartmentSubDomain.getBoundaryConditionXm().isNEUMANN()) {
Expression origExp = pdeEquation.getBoundaryXm();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryXm(new Expression(origExp + "-" + substitutedExp + "-" + diffusionRateName + "*" + exactDxName));
} else {
pdeEquation.setBoundaryXm(new Expression("-" + diffusionRateName + "*" + exactDxName));
}
} else {
throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionXm());
}
if (compartmentSubDomain.getBoundaryConditionXp().isDIRICHLET()) {
Expression origExp = pdeEquation.getBoundaryXp();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryXp(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
} else {
pdeEquation.setBoundaryXp(exactExp);
}
} else if (compartmentSubDomain.getBoundaryConditionXp().isNEUMANN()) {
Expression origExp = pdeEquation.getBoundaryXp();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryXp(new Expression(origExp + "-" + substitutedExp + "+" + diffusionRateName + "*" + exactDxName));
} else {
pdeEquation.setBoundaryXp(new Expression(diffusionRateName + "*" + exactDxName));
}
} else {
throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionXp());
}
if (compartmentSubDomain.getBoundaryConditionYm().isDIRICHLET()) {
Expression origExp = pdeEquation.getBoundaryYm();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryYm(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
} else {
pdeEquation.setBoundaryYm(exactExp);
}
} else if (compartmentSubDomain.getBoundaryConditionYm().isNEUMANN()) {
Expression origExp = pdeEquation.getBoundaryYm();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryYm(new Expression(origExp + "-" + substitutedExp + "-" + diffusionRateName + "*" + exactDyName));
} else {
pdeEquation.setBoundaryYm(new Expression("-" + diffusionRateName + "*" + exactDyName));
}
} else {
throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionYm());
}
if (compartmentSubDomain.getBoundaryConditionYp().isDIRICHLET()) {
Expression origExp = pdeEquation.getBoundaryYp();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryYp(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
} else {
pdeEquation.setBoundaryYp(exactExp);
}
} else if (compartmentSubDomain.getBoundaryConditionYp().isNEUMANN()) {
Expression origExp = pdeEquation.getBoundaryYp();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryYp(new Expression(origExp + "-" + substitutedExp + "+" + diffusionRateName + "*" + exactDyName));
} else {
pdeEquation.setBoundaryYp(new Expression(diffusionRateName + "*" + exactDyName));
}
} else {
throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionYp());
}
if (compartmentSubDomain.getBoundaryConditionZm().isDIRICHLET()) {
Expression origExp = pdeEquation.getBoundaryZm();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryZm(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
} else {
pdeEquation.setBoundaryZm(exactExp);
}
} else if (compartmentSubDomain.getBoundaryConditionZm().isNEUMANN()) {
Expression origExp = pdeEquation.getBoundaryZm();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryZm(new Expression(origExp + "-" + substitutedExp + "-" + diffusionRateName + "*" + exactDzName));
} else {
pdeEquation.setBoundaryZm(new Expression("-" + diffusionRateName + "*" + exactDzName));
}
} else {
throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionXm());
}
if (compartmentSubDomain.getBoundaryConditionZp().isDIRICHLET()) {
Expression origExp = pdeEquation.getBoundaryZp();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryZp(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
} else {
pdeEquation.setBoundaryZp(exactExp);
}
} else if (compartmentSubDomain.getBoundaryConditionZp().isNEUMANN()) {
Expression origExp = pdeEquation.getBoundaryZp();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryZp(new Expression(origExp + "-" + substitutedExp + "+" + diffusionRateName + "*" + exactDzName));
} else {
pdeEquation.setBoundaryZp(new Expression(diffusionRateName + "*" + exactDzName));
}
} else {
throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionZp());
}
} else {
throw new RuntimeException("SolverTest.constructedExactMath(): equation type " + equation.getClass().getName() + " not yet implemented");
}
}
if (subDomain instanceof MembraneSubDomain) {
MembraneSubDomain membraneSubDomain = (MembraneSubDomain) subDomain;
Enumeration<JumpCondition> enumJumpConditions = membraneSubDomain.getJumpConditions();
while (enumJumpConditions.hasMoreElements()) {
JumpCondition jumpCondition = enumJumpConditions.nextElement();
Expression origInfluxExp = jumpCondition.getInFluxExpression();
Expression origOutfluxExp = jumpCondition.getOutFluxExpression();
Expression substitutedInfluxExp = substituteWithExactSolution(origInfluxExp, membraneSubDomain, exactMath);
Expression substitutedOutfluxExp = substituteWithExactSolution(origOutfluxExp, membraneSubDomain, exactMath);
String varName = jumpCondition.getVariable().getName();
String origInfluxName = "_" + varName + "_" + subDomain.getName() + "_origInflux";
String origOutfluxName = "_" + varName + "_" + subDomain.getName() + "_origOutflux";
String substitutedInfluxName = "_" + varName + "_" + subDomain.getName() + "_substitutedInflux";
String substitutedOutfluxName = "_" + varName + "_" + subDomain.getName() + "_substitutedOutflux";
String diffusionRateInsideName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_diffusionRate";
String diffusionRateOutsideName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_diffusionRate";
String exactInsideDxName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_exact_dx";
String exactInsideDyName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_exact_dy";
String exactInsideDzName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_exact_dz";
String exactOutsideDxName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_exact_dx";
String exactOutsideDyName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_exact_dy";
String exactOutsideDzName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_exact_dz";
String outwardNormalXName = "_" + membraneSubDomain.getInsideCompartment().getName() + "_Nx";
String outwardNormalYName = "_" + membraneSubDomain.getInsideCompartment().getName() + "_Ny";
String outwardNormalZName = "_" + membraneSubDomain.getInsideCompartment().getName() + "_Nz";
String exactInfluxName = "_" + varName + "_" + membraneSubDomain.getName() + "_exactInflux";
String exactOutfluxName = "_" + varName + "_" + membraneSubDomain.getName() + "_exactOutflux";
Expression exactInfluxExp = new Expression(diffusionRateInsideName + " * (" + outwardNormalXName + "*" + exactInsideDxName + " + " + outwardNormalYName + "*" + exactInsideDyName + " + " + outwardNormalZName + "*" + exactInsideDzName + ")");
Expression exactOutfluxExp = new Expression("-" + diffusionRateOutsideName + " * (" + outwardNormalXName + "*" + exactOutsideDxName + " + " + outwardNormalYName + "*" + exactOutsideDyName + " + " + outwardNormalZName + "*" + exactOutsideDzName + ")");
Expression newInfluxExp = new Expression(origInfluxName + " - " + substitutedInfluxName + " + " + exactInfluxName);
Expression newOutfluxExp = new Expression(origOutfluxName + " - " + substitutedOutfluxName + " + " + exactOutfluxName);
varHash.addVariable(new Function(origInfluxName, origInfluxExp, domain));
varHash.addVariable(new Function(origOutfluxName, origOutfluxExp, domain));
varHash.addVariable(new Function(exactInfluxName, exactInfluxExp, domain));
varHash.addVariable(new Function(exactOutfluxName, exactOutfluxExp, domain));
varHash.addVariable(new Function(substitutedInfluxName, substitutedInfluxExp, domain));
varHash.addVariable(new Function(substitutedOutfluxName, substitutedOutfluxExp, domain));
jumpCondition.setInFlux(newInfluxExp);
jumpCondition.setOutFlux(newOutfluxExp);
}
}
}
exactMath.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
if (!exactMath.isValid()) {
throw new RuntimeException("generated Math is not valid: " + exactMath.getWarning());
}
return exactMath;
}
Also used :
JumpCondition(cbit.vcell.math.JumpCondition)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable)
InsideVariable(cbit.vcell.math.InsideVariable)
SensVariable(cbit.vcell.solver.ode.SensVariable)
FilamentVariable(cbit.vcell.math.FilamentVariable)
VolVariable(cbit.vcell.math.VolVariable)
VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable)
ReservedVariable(cbit.vcell.math.ReservedVariable)
MemVariable(cbit.vcell.math.MemVariable)
OutsideVariable(cbit.vcell.math.OutsideVariable)
FilamentRegionVariable(cbit.vcell.math.FilamentRegionVariable)
Variable(cbit.vcell.math.Variable)
MathDescription(cbit.vcell.math.MathDescription)
VariableHash(cbit.vcell.math.VariableHash)
Constant(cbit.vcell.math.Constant)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
PdeEquation(cbit.vcell.math.PdeEquation)
Function(cbit.vcell.math.Function)
OdeEquation(cbit.vcell.math.OdeEquation)
PdeEquation(cbit.vcell.math.PdeEquation)
Equation(cbit.vcell.math.Equation)
SolutionTemplate(cbit.vcell.numericstest.SolutionTemplate)
ConstructedSolutionTemplate(cbit.vcell.numericstest.ConstructedSolutionTemplate)
OdeEquation(cbit.vcell.math.OdeEquation)
Expression(cbit.vcell.parser.Expression)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
Domain(cbit.vcell.math.Variable.Domain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume)
Aggregations
MembraneSubDomain (cbit.vcell.math.MembraneSubDomain)31
CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)26
SubDomain (cbit.vcell.math.SubDomain)21
Expression (cbit.vcell.parser.Expression)13
ArrayList (java.util.ArrayList)13
MathDescription (cbit.vcell.math.MathDescription)12
SubVolume (cbit.vcell.geometry.SubVolume)10
ExpressionException (cbit.vcell.parser.ExpressionException)10
Variable (cbit.vcell.math.Variable)9
SurfaceClass (cbit.vcell.geometry.SurfaceClass)8
Constant (cbit.vcell.math.Constant)7
Equation (cbit.vcell.math.Equation)7
MathException (cbit.vcell.math.MathException)7
OdeEquation (cbit.vcell.math.OdeEquation)7
JumpCondition (cbit.vcell.math.JumpCondition)6
ParticleProperties (cbit.vcell.math.ParticleProperties)6
HashMap (java.util.HashMap)6
AnalyticSubVolume (cbit.vcell.geometry.AnalyticSubVolume)5
Function (cbit.vcell.math.Function)5
MemVariable (cbit.vcell.math.MemVariable)5
