Example 16 with SurfaceClass
use of cbit.vcell.geometry.SurfaceClass in project vcell by virtualcell.
the class DiffEquMathMapping method addSpatialProcesses.
private void addSpatialProcesses(VariableHash varHash, ArrayList<CompartmentSubdomainContext> compartmentSubdomainContexts, ArrayList<MembraneSubdomainContext> membraneSubdomainContexts) throws MathException, MappingException, ExpressionException {
if (simContext.getGeometry().getDimension() == 0) {
return;
}
//
for (SpatialObject spatialObject : simContext.getSpatialObjects()) {
if (spatialObject instanceof PointObject) {
PointObject pointObject = (PointObject) spatialObject;
//
// if true, have to solve for this category
//
boolean bPosition = pointObject.isQuantityCategoryEnabled(QuantityCategory.PointPosition);
boolean bVelocity = pointObject.isQuantityCategoryEnabled(QuantityCategory.PointVelocity);
boolean bDirection = pointObject.isQuantityCategoryEnabled(QuantityCategory.DirectionToPoint);
boolean bDistance = pointObject.isQuantityCategoryEnabled(QuantityCategory.PointDistanceMap);
//
// either make a point subdomain, or just define functions.
//
ArrayList<PointLocation> pointLocationProcesses = new ArrayList<PointLocation>();
ArrayList<PointKinematics> pointKinematicsProcesses = new ArrayList<PointKinematics>();
for (SpatialProcess spatialProcess : simContext.getSpatialProcesses()) {
if (spatialProcess instanceof PointLocation && ((PointLocation) spatialProcess).getPointObject() == pointObject) {
pointLocationProcesses.add((PointLocation) spatialProcess);
}
if (spatialProcess instanceof PointKinematics && ((PointKinematics) spatialProcess).getPointObject() == pointObject) {
pointKinematicsProcesses.add((PointKinematics) spatialProcess);
}
}
if (pointLocationProcesses.size() == 1 && pointKinematicsProcesses.size() == 0 && !bVelocity) {
GeometryClass gc = null;
PointLocation pointLocation = pointLocationProcesses.get(0);
if (bPosition) {
if (simContext.getGeometry().getDimension() == 1) {
SpatialQuantity posXQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.X);
LocalParameter posXParam = pointLocation.getParameter(SpatialProcessParameterType.PointPositionX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posXParam, gc), getIdentifierSubstitutions(posXParam.getExpression(), posXParam.getUnitDefinition(), gc), gc));
Expression posXExp = new Expression(posXParam, pointLocation.getNameScope());
Expression xExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.X), simContext.getModel().getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posXQuantity, gc), getIdentifierSubstitutions(posXExp, posXQuantity.getUnitDefinition(), gc), gc));
Expression posX_minus_X = Expression.add(posXExp, Expression.negate(xExp));
Expression signum_posX_minus_X = Expression.function(FunctionType.MIN, new Expression(1.0), Expression.function(FunctionType.MAX, Expression.mult(new Expression(1e10), posX_minus_X), new Expression(-1)));
if (bDirection) {
SpatialQuantity dirXQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.X);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirXQuantity, gc), getIdentifierSubstitutions(signum_posX_minus_X, dirXQuantity.getUnitDefinition(), gc), gc));
}
if (bDistance) {
Expression abs_X_minux_posX = Expression.function(FunctionType.ABS, posX_minus_X);
SpatialQuantity distanceQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointDistanceMap, QuantityComponent.Scalar);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(distanceQuantity, gc), getIdentifierSubstitutions(abs_X_minux_posX, distanceQuantity.getUnitDefinition(), gc), gc));
}
}
if (simContext.getGeometry().getDimension() == 2) {
SpatialQuantity posXQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.X);
LocalParameter posXParam = pointLocation.getParameter(SpatialProcessParameterType.PointPositionX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posXParam, gc), getIdentifierSubstitutions(posXParam.getExpression(), posXParam.getUnitDefinition(), gc), gc));
Expression posXExp = new Expression(posXParam, pointLocation.getNameScope());
Expression xExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.X), simContext.getModel().getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posXQuantity, gc), getIdentifierSubstitutions(new Expression(posXParam, pointLocation.getNameScope()), posXQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity posYQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.Y);
LocalParameter posYParam = pointLocation.getParameter(SpatialProcessParameterType.PointPositionY);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posYParam, gc), getIdentifierSubstitutions(posYParam.getExpression(), posYParam.getUnitDefinition(), gc), gc));
Expression posYExp = new Expression(posYParam, pointLocation.getNameScope());
Expression yExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.Y), simContext.getModel().getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posYQuantity, gc), getIdentifierSubstitutions(new Expression(posYParam, pointLocation.getNameScope()), posYQuantity.getUnitDefinition(), gc), gc));
Expression posX_minux_X = Expression.add(posXExp, Expression.negate(xExp));
Expression posY_minux_Y = Expression.add(posYExp, Expression.negate(yExp));
Expression DX2 = Expression.mult(posX_minux_X, posX_minux_X);
Expression DY2 = Expression.mult(posY_minux_Y, posY_minux_Y);
Expression sqrt_DX2_DY2 = Expression.function(FunctionType.SQRT, Expression.add(DX2, DY2));
Expression dirX = Expression.div(posX_minux_X, Expression.add(sqrt_DX2_DY2, new Expression(1e-8)));
Expression dirY = Expression.div(posY_minux_Y, Expression.add(sqrt_DX2_DY2, new Expression(1e-8)));
if (bDirection) {
SpatialQuantity dirXQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.X);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirXQuantity, gc), getIdentifierSubstitutions(dirX, dirXQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity dirYQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.Y);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirYQuantity, gc), getIdentifierSubstitutions(dirY, dirYQuantity.getUnitDefinition(), gc), gc));
}
if (bDistance) {
SpatialQuantity distanceQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointDistanceMap, QuantityComponent.Scalar);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(distanceQuantity, gc), getIdentifierSubstitutions(sqrt_DX2_DY2, distanceQuantity.getUnitDefinition(), gc), gc));
}
}
if (simContext.getGeometry().getDimension() == 3) {
SpatialQuantity posXQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.X);
LocalParameter posXParam = pointLocation.getParameter(SpatialProcessParameterType.PointPositionX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posXParam, gc), getIdentifierSubstitutions(posXParam.getExpression(), posXParam.getUnitDefinition(), gc), gc));
Expression posXExp = new Expression(posXParam, pointLocation.getNameScope());
Expression xExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.X), simContext.getModel().getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posXQuantity, gc), getIdentifierSubstitutions(new Expression(posXParam, pointLocation.getNameScope()), posXQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity posYQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.Y);
LocalParameter posYParam = pointLocation.getParameter(SpatialProcessParameterType.PointPositionY);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posYParam, gc), getIdentifierSubstitutions(posYParam.getExpression(), posYParam.getUnitDefinition(), gc), gc));
Expression posYExp = new Expression(posYParam, pointLocation.getNameScope());
Expression yExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.Y), simContext.getModel().getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posYQuantity, gc), getIdentifierSubstitutions(new Expression(posYParam, pointLocation.getNameScope()), posYQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity posZQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.Z);
LocalParameter posZParam = pointLocation.getParameter(SpatialProcessParameterType.PointPositionZ);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posZParam, gc), getIdentifierSubstitutions(posZParam.getExpression(), posZParam.getUnitDefinition(), gc), gc));
Expression posZExp = new Expression(posZParam, pointLocation.getNameScope());
Expression zExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.Z), simContext.getModel().getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posZQuantity, gc), getIdentifierSubstitutions(new Expression(posZParam, pointLocation.getNameScope()), posZQuantity.getUnitDefinition(), gc), gc));
Expression posX_minux_X = Expression.add(posXExp, Expression.negate(xExp));
Expression posY_minux_Y = Expression.add(posYExp, Expression.negate(yExp));
Expression posZ_minux_Z = Expression.add(posZExp, Expression.negate(zExp));
Expression DX2 = Expression.mult(posX_minux_X, posX_minux_X);
Expression DY2 = Expression.mult(posY_minux_Y, posY_minux_Y);
Expression DZ2 = Expression.mult(posZ_minux_Z, posZ_minux_Z);
Expression sqrt_DX2_DY2_DZ2 = Expression.function(FunctionType.SQRT, Expression.add(DX2, DY2, DZ2));
Expression dirX = Expression.div(posX_minux_X, Expression.add(sqrt_DX2_DY2_DZ2, new Expression(1e-8)));
Expression dirY = Expression.div(posY_minux_Y, Expression.add(sqrt_DX2_DY2_DZ2, new Expression(1e-8)));
Expression dirZ = Expression.div(posZ_minux_Z, Expression.add(sqrt_DX2_DY2_DZ2, new Expression(1e-8)));
if (bDirection) {
SpatialQuantity dirXQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.X);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirXQuantity, gc), getIdentifierSubstitutions(dirX, dirXQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity dirYQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.Y);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirYQuantity, gc), getIdentifierSubstitutions(dirY, dirYQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity dirZQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.Z);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirZQuantity, gc), getIdentifierSubstitutions(dirZ, dirZQuantity.getUnitDefinition(), gc), gc));
}
if (bDistance) {
SpatialQuantity distanceQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointDistanceMap, QuantityComponent.Scalar);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(distanceQuantity, gc), getIdentifierSubstitutions(sqrt_DX2_DY2_DZ2, distanceQuantity.getUnitDefinition(), gc), gc));
}
}
} else {
throw new MappingException("PointLocation process defined for pointObject '" + pointObject.getName() + "' but Position not enabled");
}
} else if (pointLocationProcesses.size() == 0 && pointKinematicsProcesses.size() == 1) {
GeometryClass gc = null;
PointKinematics pointKinematics = pointKinematicsProcesses.get(0);
if (bPosition && bVelocity) {
LocalParameter velXParam = pointKinematics.getParameter(SpatialProcessParameterType.PointVelocityX);
LocalParameter iniPosXParam = pointKinematics.getParameter(SpatialProcessParameterType.PointInitialPositionX);
SpatialQuantity posXQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.X);
SpatialQuantity velXQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointVelocity, QuantityComponent.X);
Expression posXExp = new Expression(posXQuantity, simContext.getNameScope());
Expression xExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.X), simContext.getModel().getNameScope());
Expression posX_minux_X = Expression.add(posXExp, Expression.negate(xExp));
LocalParameter velYParam = pointKinematics.getParameter(SpatialProcessParameterType.PointVelocityY);
LocalParameter iniPosYParam = pointKinematics.getParameter(SpatialProcessParameterType.PointInitialPositionY);
SpatialQuantity posYQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.Y);
SpatialQuantity velYQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointVelocity, QuantityComponent.Y);
Expression posYExp = new Expression(posYQuantity, simContext.getNameScope());
Expression yExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.Y), simContext.getModel().getNameScope());
Expression posY_minux_Y = Expression.add(posYExp, Expression.negate(yExp));
LocalParameter velZParam = pointKinematics.getParameter(SpatialProcessParameterType.PointVelocityZ);
LocalParameter iniPosZParam = pointKinematics.getParameter(SpatialProcessParameterType.PointInitialPositionZ);
SpatialQuantity posZQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.Z);
SpatialQuantity velZQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointVelocity, QuantityComponent.Z);
Expression posZExp = new Expression(posZQuantity, simContext.getNameScope());
Expression zExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.Z), simContext.getModel().getNameScope());
Expression posZ_minux_Z = Expression.add(posZExp, Expression.negate(zExp));
String pointSubdomainName = pointObject.getName();
Domain domain = new Domain(pointSubdomainName);
PointSubDomain pointSubdomain = new PointSubDomain(pointSubdomainName);
mathDesc.addSubDomain(pointSubdomain);
if (simContext.getGeometry().getDimension() >= 1) {
PointVariable posXVar = new PointVariable(getMathSymbol(posXQuantity, gc), domain);
varHash.addVariable(posXVar);
Expression initXExp = getIdentifierSubstitutions(new Expression(iniPosXParam, pointKinematics.getNameScope()), iniPosXParam.getUnitDefinition(), gc);
Expression rateXExp = getIdentifierSubstitutions(new Expression(velXParam, pointKinematics.getNameScope()), velXParam.getUnitDefinition(), gc);
OdeEquation odeX = new OdeEquation(posXVar, initXExp, rateXExp);
pointSubdomain.addEquation(odeX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(iniPosXParam, gc), getIdentifierSubstitutions(new Expression(iniPosXParam.getExpression()), iniPosXParam.getUnitDefinition(), gc), gc));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velXParam, gc), getIdentifierSubstitutions(new Expression(velXParam.getExpression()), velXParam.getUnitDefinition(), gc), gc));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velXQuantity, gc), getIdentifierSubstitutions(new Expression(velXParam, pointKinematics.getNameScope()), velXQuantity.getUnitDefinition(), gc), gc));
pointSubdomain.setPositionX(getIdentifierSubstitutions(new Expression(posXQuantity, simContext.getNameScope()), posXQuantity.getUnitDefinition(), gc));
}
if (simContext.getGeometry().getDimension() >= 2) {
PointVariable posYVar = new PointVariable(getMathSymbol(posYQuantity, gc), domain);
varHash.addVariable(posYVar);
Expression initYExp = getIdentifierSubstitutions(new Expression(iniPosYParam, pointKinematics.getNameScope()), iniPosYParam.getUnitDefinition(), gc);
Expression rateYExp = getIdentifierSubstitutions(new Expression(velYParam, pointKinematics.getNameScope()), velYParam.getUnitDefinition(), gc);
OdeEquation odeY = new OdeEquation(posYVar, initYExp, rateYExp);
pointSubdomain.addEquation(odeY);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(iniPosYParam, gc), getIdentifierSubstitutions(new Expression(iniPosYParam.getExpression()), iniPosYParam.getUnitDefinition(), gc), gc));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velYParam, gc), getIdentifierSubstitutions(new Expression(velYParam.getExpression()), velYParam.getUnitDefinition(), gc), gc));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velYQuantity, gc), getIdentifierSubstitutions(new Expression(velYParam, pointKinematics.getNameScope()), velYQuantity.getUnitDefinition(), gc), gc));
pointSubdomain.setPositionY(getIdentifierSubstitutions(new Expression(posYQuantity, simContext.getNameScope()), posYQuantity.getUnitDefinition(), gc));
}
if (simContext.getGeometry().getDimension() == 3) {
PointVariable posZVar = new PointVariable(getMathSymbol(posZQuantity, gc), domain);
varHash.addVariable(posZVar);
Expression initZExp = getIdentifierSubstitutions(new Expression(iniPosZParam, pointKinematics.getNameScope()), iniPosZParam.getUnitDefinition(), gc);
Expression rateZExp = getIdentifierSubstitutions(new Expression(velZParam, pointKinematics.getNameScope()), velZParam.getUnitDefinition(), gc);
OdeEquation odeZ = new OdeEquation(posZVar, initZExp, rateZExp);
pointSubdomain.addEquation(odeZ);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(iniPosZParam, gc), getIdentifierSubstitutions(new Expression(iniPosZParam.getExpression()), iniPosZParam.getUnitDefinition(), gc), gc));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velZParam, gc), getIdentifierSubstitutions(new Expression(velZParam.getExpression()), velZParam.getUnitDefinition(), gc), gc));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velZQuantity, gc), getIdentifierSubstitutions(new Expression(velZParam, pointKinematics.getNameScope()), velZQuantity.getUnitDefinition(), gc), gc));
pointSubdomain.setPositionZ(getIdentifierSubstitutions(new Expression(posZQuantity, simContext.getNameScope()), posZQuantity.getUnitDefinition(), gc));
}
if (simContext.getGeometry().getDimension() == 1) {
Expression signum_posX_minus_X = Expression.function(FunctionType.MIN, new Expression(1.0), Expression.function(FunctionType.MAX, Expression.mult(new Expression(1e10), posX_minux_X), new Expression(-1)));
if (bDirection) {
SpatialQuantity dirXQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.X);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirXQuantity, gc), getIdentifierSubstitutions(signum_posX_minus_X, dirXQuantity.getUnitDefinition(), gc), gc));
}
if (bDistance) {
Expression abs_X_minux_posX = Expression.function(FunctionType.ABS, posX_minux_X);
SpatialQuantity distanceQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointDistanceMap, QuantityComponent.Scalar);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(distanceQuantity, gc), getIdentifierSubstitutions(abs_X_minux_posX, distanceQuantity.getUnitDefinition(), gc), gc));
}
}
if (simContext.getGeometry().getDimension() == 2) {
Expression DX2 = Expression.mult(posX_minux_X, posX_minux_X);
Expression DY2 = Expression.mult(posY_minux_Y, posY_minux_Y);
Expression sqrt_DX2_DY2 = Expression.function(FunctionType.SQRT, Expression.add(DX2, DY2));
Expression dirX = Expression.div(posX_minux_X, Expression.add(sqrt_DX2_DY2, new Expression(1e-8)));
Expression dirY = Expression.div(posY_minux_Y, Expression.add(sqrt_DX2_DY2, new Expression(1e-8)));
if (bDirection) {
SpatialQuantity dirXQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.X);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirXQuantity, gc), getIdentifierSubstitutions(dirX, dirXQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity dirYQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.Y);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirYQuantity, gc), getIdentifierSubstitutions(dirY, dirYQuantity.getUnitDefinition(), gc), gc));
}
if (bDistance) {
SpatialQuantity distanceQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointDistanceMap, QuantityComponent.Scalar);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(distanceQuantity, gc), getIdentifierSubstitutions(sqrt_DX2_DY2, distanceQuantity.getUnitDefinition(), gc), gc));
}
}
if (simContext.getGeometry().getDimension() == 3) {
Expression DX2 = Expression.mult(posX_minux_X, posX_minux_X);
Expression DY2 = Expression.mult(posY_minux_Y, posY_minux_Y);
Expression DZ2 = Expression.mult(posZ_minux_Z, posZ_minux_Z);
Expression sqrt_DX2_DY2_DZ2 = Expression.function(FunctionType.SQRT, Expression.add(DX2, DY2, DZ2));
Expression dirX = Expression.div(posX_minux_X, Expression.add(sqrt_DX2_DY2_DZ2, new Expression(1e-8)));
Expression dirY = Expression.div(posY_minux_Y, Expression.add(sqrt_DX2_DY2_DZ2, new Expression(1e-8)));
Expression dirZ = Expression.div(posZ_minux_Z, Expression.add(sqrt_DX2_DY2_DZ2, new Expression(1e-8)));
if (bDirection) {
SpatialQuantity dirXQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.X);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirXQuantity, gc), getIdentifierSubstitutions(dirX, dirXQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity dirYQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.Y);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirYQuantity, gc), getIdentifierSubstitutions(dirY, dirYQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity dirZQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.Z);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirZQuantity, gc), getIdentifierSubstitutions(dirZ, dirZQuantity.getUnitDefinition(), gc), gc));
}
if (bDistance) {
SpatialQuantity distanceQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointDistanceMap, QuantityComponent.Scalar);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(distanceQuantity, gc), getIdentifierSubstitutions(sqrt_DX2_DY2_DZ2, distanceQuantity.getUnitDefinition(), gc), gc));
}
}
} else {
throw new MappingException(pointKinematics.getDescription() + " process defined for pointObject '" + pointObject.getName() + "' but Position and Velocity not enabled");
}
} else {
throw new MappingException("expecting 1 location or kinematics process for point '" + pointObject.getName() + "'");
}
} else if (spatialObject instanceof SurfaceRegionObject) {
SurfaceRegionObject surfaceRegionObject = (SurfaceRegionObject) spatialObject;
SubVolume insideSubvolume = surfaceRegionObject.getInsideSubVolume();
SubVolume outsideSubvolume = surfaceRegionObject.getOutsideSubVolume();
SurfaceClass surfaceClass = simContext.getGeometry().getGeometrySurfaceDescription().getSurfaceClass(insideSubvolume, outsideSubvolume);
MembraneSubdomainContext memSubdomainContext = null;
for (MembraneSubdomainContext context : membraneSubdomainContexts) {
if (context.surfaceClass == surfaceClass) {
memSubdomainContext = context;
}
}
//
// if true, have to solve for this category
//
boolean bNormal = surfaceRegionObject.isQuantityCategoryEnabled(QuantityCategory.Normal);
boolean bVelocity = surfaceRegionObject.isQuantityCategoryEnabled(QuantityCategory.SurfaceVelocity);
boolean bDistance = surfaceRegionObject.isQuantityCategoryEnabled(QuantityCategory.SurfaceDistanceMap);
boolean bDirection = surfaceRegionObject.isQuantityCategoryEnabled(QuantityCategory.DirectionToSurface);
boolean bSize = surfaceRegionObject.isQuantityCategoryEnabled(QuantityCategory.SurfaceSize);
if (bVelocity) {
ArrayList<SurfaceKinematics> surfaceKinematicsList = new ArrayList<SurfaceKinematics>();
for (SpatialProcess spatialProcess : simContext.getSpatialProcesses()) {
if (spatialProcess instanceof SurfaceKinematics && ((SurfaceKinematics) spatialProcess).getSurfaceRegionObject() == surfaceRegionObject) {
surfaceKinematicsList.add((SurfaceKinematics) spatialProcess);
}
}
if (surfaceKinematicsList.size() == 1) {
SurfaceKinematics surfaceKinematics = surfaceKinematicsList.get(0);
if (simContext.getGeometry().getDimension() >= 1) {
SpatialQuantity velXQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.X);
LocalParameter velXParam = surfaceKinematics.getParameter(SpatialProcessParameterType.SurfaceVelocityX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velXParam, surfaceClass), getIdentifierSubstitutions(velXParam.getExpression(), velXParam.getUnitDefinition(), surfaceClass), surfaceClass));
Expression velXExp = new Expression(velXParam, surfaceKinematics.getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velXQuantity, surfaceClass), getIdentifierSubstitutions(velXExp, velXQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
if (bNormal) {
SpatialQuantity normXQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.Normal, QuantityComponent.X);
Expression normXExp = new Expression(MathFunctionDefinitions.FUNCTION_normalX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(normXQuantity, surfaceClass), getIdentifierSubstitutions(normXExp, normXQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
}
}
if (simContext.getGeometry().getDimension() >= 2) {
SpatialQuantity velYQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.Y);
LocalParameter velYParam = surfaceKinematics.getParameter(SpatialProcessParameterType.SurfaceVelocityY);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velYParam, surfaceClass), getIdentifierSubstitutions(velYParam.getExpression(), velYParam.getUnitDefinition(), surfaceClass), surfaceClass));
Expression velYExp = new Expression(velYParam, surfaceKinematics.getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velYQuantity, surfaceClass), getIdentifierSubstitutions(velYExp, velYQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
if (bNormal) {
SpatialQuantity normYQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.Normal, QuantityComponent.Y);
Expression normYExp = new Expression(MathFunctionDefinitions.FUNCTION_normalY);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(normYQuantity, surfaceClass), getIdentifierSubstitutions(normYExp, normYQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
}
}
if (simContext.getGeometry().getDimension() >= 3) {
SpatialQuantity velZQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.Z);
LocalParameter velZParam = surfaceKinematics.getParameter(SpatialProcessParameterType.SurfaceVelocityZ);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velZParam, surfaceClass), getIdentifierSubstitutions(velZParam.getExpression(), velZParam.getUnitDefinition(), surfaceClass), surfaceClass));
Expression velYExp = new Expression(velZParam, surfaceKinematics.getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velZQuantity, surfaceClass), getIdentifierSubstitutions(velYExp, velZQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
if (bNormal) {
SpatialQuantity normZQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.Normal, QuantityComponent.Z);
// MathFunctionDefinitions.FUNCTION_normalZ);
Expression normZExp = new Expression("normalZ_not_implemented()");
varHash.addVariable(newFunctionOrConstant(getMathSymbol(normZQuantity, surfaceClass), getIdentifierSubstitutions(normZExp, normZQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
}
}
} else {
throw new MappingException("expecting 1 Surface Kinematics process for Surface Object '" + surfaceRegionObject.getName() + "'");
}
}
if (bSize) {
SpatialQuantity sizeQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceSize, QuantityComponent.Scalar);
String funcName = MathFunctionDefinitions.Function_regionArea_current.getFunctionName();
Expression sizeExp = Expression.function(funcName, new Expression[] { new Expression("'" + surfaceClass.getName() + "'") });
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeQuantity, surfaceClass), getIdentifierSubstitutions(sizeExp, sizeQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
}
if (bDirection) {
SpatialQuantity directionXQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.DirectionToSurface, QuantityComponent.X);
SpatialQuantity directionYQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.DirectionToSurface, QuantityComponent.Y);
SpatialQuantity directionZQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.DirectionToSurface, QuantityComponent.Z);
for (SubVolume adjacentSubvolume : surfaceClass.getAdjacentSubvolumes()) {
Domain domain = new Domain(adjacentSubvolume);
CompartmentSubDomain compSubdomain = mathDesc.getCompartmentSubDomain(adjacentSubvolume.getName());
if (simContext.getGeometry().getDimension() >= 1) {
String name = adjacentSubvolume.getName() + "_dirX_" + memSubdomainContext.membraneSubdomain.getName();
LocalizedDirectionToMembraneQuantity localDirectionXQuantity = addLocalizedDirectionToMembraneQuantity(name, surfaceClass, adjacentSubvolume, QuantityComponent.X);
VolVariable distanceVar = new VolVariable(getMathSymbol(localDirectionXQuantity, adjacentSubvolume), domain);
varHash.addVariable(distanceVar);
ComputeMembraneMetricEquation membraneMetricEquation = new ComputeMembraneMetricEquation(distanceVar, MembraneMetricComponent.directionToMembraneX);
membraneMetricEquation.setTargetMembraneName(memSubdomainContext.membraneSubdomain.getName());
compSubdomain.addEquation(membraneMetricEquation);
}
if (simContext.getGeometry().getDimension() >= 2) {
String name = adjacentSubvolume.getName() + "_dirY_" + memSubdomainContext.membraneSubdomain.getName();
LocalizedDirectionToMembraneQuantity localDirectionYQuantity = addLocalizedDirectionToMembraneQuantity(name, surfaceClass, adjacentSubvolume, QuantityComponent.Y);
VolVariable distanceVar = new VolVariable(getMathSymbol(localDirectionYQuantity, adjacentSubvolume), domain);
varHash.addVariable(distanceVar);
ComputeMembraneMetricEquation membraneMetricEquation = new ComputeMembraneMetricEquation(distanceVar, MembraneMetricComponent.directionToMembraneY);
membraneMetricEquation.setTargetMembraneName(memSubdomainContext.membraneSubdomain.getName());
compSubdomain.addEquation(membraneMetricEquation);
}
if (simContext.getGeometry().getDimension() == 3) {
String name = adjacentSubvolume.getName() + "_dirZ_" + memSubdomainContext.membraneSubdomain.getName();
LocalizedDirectionToMembraneQuantity localDirectionZQuantity = addLocalizedDirectionToMembraneQuantity(name, surfaceClass, adjacentSubvolume, QuantityComponent.Z);
VolVariable distanceVar = new VolVariable(getMathSymbol(localDirectionZQuantity, adjacentSubvolume), domain);
varHash.addVariable(distanceVar);
ComputeMembraneMetricEquation membraneMetricEquation = new ComputeMembraneMetricEquation(distanceVar, MembraneMetricComponent.directionToMembraneZ);
membraneMetricEquation.setTargetMembraneName(memSubdomainContext.membraneSubdomain.getName());
compSubdomain.addEquation(membraneMetricEquation);
}
}
}
if (bDistance) {
SpatialQuantity distanceQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceDistanceMap, QuantityComponent.Scalar);
for (SubVolume adjacentSubvolume : surfaceClass.getAdjacentSubvolumes()) {
String name = adjacentSubvolume.getName() + "_distanceTo_" + memSubdomainContext.membraneSubdomain.getName();
LocalizedDistanceToMembraneQuantity localDistanceQuantity = addLocalizedDistanceToMembraneQuantity(name, surfaceClass, adjacentSubvolume);
CompartmentSubDomain compSubdomain = mathDesc.getCompartmentSubDomain(adjacentSubvolume.getName());
Domain domain = new Domain(adjacentSubvolume);
VolVariable distanceVar = new VolVariable(getMathSymbol(localDistanceQuantity, adjacentSubvolume), domain);
varHash.addVariable(distanceVar);
ComputeMembraneMetricEquation membraneMetricEquation = new ComputeMembraneMetricEquation(distanceVar, MembraneMetricComponent.distanceToMembrane);
membraneMetricEquation.setTargetMembraneName(memSubdomainContext.membraneSubdomain.getName());
compSubdomain.addEquation(membraneMetricEquation);
}
}
} else if (spatialObject instanceof VolumeRegionObject) {
VolumeRegionObject volumeRegionObject = (VolumeRegionObject) spatialObject;
SubVolume subvolume = volumeRegionObject.getSubVolume();
boolean bCentroid = volumeRegionObject.isQuantityCategoryEnabled(QuantityCategory.DirectionToSurface);
boolean bSize = volumeRegionObject.isQuantityCategoryEnabled(QuantityCategory.VolumeSize);
boolean bVelocity = volumeRegionObject.isQuantityCategoryEnabled(QuantityCategory.InteriorVelocity);
if (bVelocity) {
ArrayList<VolumeKinematics> volumeKinematicsList = new ArrayList<VolumeKinematics>();
for (SpatialProcess spatialProcess : simContext.getSpatialProcesses()) {
if (spatialProcess instanceof VolumeKinematics && ((VolumeKinematics) spatialProcess).getVolumeRegionObject() == volumeRegionObject) {
volumeKinematicsList.add((VolumeKinematics) spatialProcess);
}
}
if (volumeKinematicsList.size() == 1) {
VolumeKinematics volumeKinematics = volumeKinematicsList.get(0);
if (simContext.getGeometry().getDimension() >= 1) {
SpatialQuantity velXQuantity = volumeRegionObject.getSpatialQuantity(QuantityCategory.InteriorVelocity, QuantityComponent.X);
LocalParameter velXParam = volumeKinematics.getParameter(SpatialProcessParameterType.InternalVelocityX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velXParam, subvolume), getIdentifierSubstitutions(velXParam.getExpression(), velXParam.getUnitDefinition(), subvolume), subvolume));
Expression velXExp = new Expression(velXParam, volumeKinematics.getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velXQuantity, subvolume), getIdentifierSubstitutions(velXExp, velXQuantity.getUnitDefinition(), subvolume), subvolume));
}
if (simContext.getGeometry().getDimension() >= 2) {
SpatialQuantity velYQuantity = volumeRegionObject.getSpatialQuantity(QuantityCategory.InteriorVelocity, QuantityComponent.Y);
LocalParameter velYParam = volumeKinematics.getParameter(SpatialProcessParameterType.InternalVelocityY);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velYParam, subvolume), getIdentifierSubstitutions(velYParam.getExpression(), velYParam.getUnitDefinition(), subvolume), subvolume));
Expression velYExp = new Expression(velYParam, volumeKinematics.getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velYQuantity, subvolume), getIdentifierSubstitutions(velYExp, velYQuantity.getUnitDefinition(), subvolume), subvolume));
}
if (simContext.getGeometry().getDimension() == 3) {
SpatialQuantity velZQuantity = volumeRegionObject.getSpatialQuantity(QuantityCategory.InteriorVelocity, QuantityComponent.Z);
LocalParameter velZParam = volumeKinematics.getParameter(SpatialProcessParameterType.InternalVelocityZ);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velZParam, subvolume), getIdentifierSubstitutions(velZParam.getExpression(), velZParam.getUnitDefinition(), subvolume), subvolume));
Expression velZExp = new Expression(velZParam, volumeKinematics.getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velZQuantity, subvolume), getIdentifierSubstitutions(velZExp, velZQuantity.getUnitDefinition(), subvolume), subvolume));
}
} else {
throw new MappingException("expecting 1 Volume Kinematics process for Volume Object '" + volumeRegionObject.getName() + "'");
}
}
if (bSize) {
SpatialQuantity sizeQuantity = volumeRegionObject.getSpatialQuantity(QuantityCategory.VolumeSize, QuantityComponent.Scalar);
String funcName = MathFunctionDefinitions.Function_regionVolume_current.getFunctionName();
Expression sizeExp = Expression.function(funcName, new Expression[] { new Expression("'" + subvolume.getName() + "'") });
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeQuantity, subvolume), getIdentifierSubstitutions(sizeExp, sizeQuantity.getUnitDefinition(), subvolume), subvolume));
}
if (bCentroid) {
throw new MappingException(QuantityCategory.Centroid.description + " not yet implemented in math generation");
}
}
}
}
Also used :
GeometryClass(cbit.vcell.geometry.GeometryClass)
PointLocation(cbit.vcell.mapping.spatial.processes.PointLocation)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
VolumeKinematics(cbit.vcell.mapping.spatial.processes.VolumeKinematics)
ArrayList(java.util.ArrayList)
PointSubDomain(cbit.vcell.math.PointSubDomain)
SpatialObject(cbit.vcell.mapping.spatial.SpatialObject)
PointObject(cbit.vcell.mapping.spatial.PointObject)
SpatialProcess(cbit.vcell.mapping.spatial.processes.SpatialProcess)
SubVolume(cbit.vcell.geometry.SubVolume)
ComputeMembraneMetricEquation(cbit.vcell.math.ComputeMembraneMetricEquation)
SpatialQuantity(cbit.vcell.mapping.spatial.SpatialObject.SpatialQuantity)
VolumeRegionObject(cbit.vcell.mapping.spatial.VolumeRegionObject)
VolVariable(cbit.vcell.math.VolVariable)
SurfaceKinematics(cbit.vcell.mapping.spatial.processes.SurfaceKinematics)
LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter)
OdeEquation(cbit.vcell.math.OdeEquation)
Expression(cbit.vcell.parser.Expression)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
PointKinematics(cbit.vcell.mapping.spatial.processes.PointKinematics)
PointVariable(cbit.vcell.math.PointVariable)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
PointSubDomain(cbit.vcell.math.PointSubDomain)
Domain(cbit.vcell.math.Variable.Domain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
SurfaceRegionObject(cbit.vcell.mapping.spatial.SurfaceRegionObject)
Example 17 with SurfaceClass
use of cbit.vcell.geometry.SurfaceClass in project vcell by virtualcell.
the class FeatureMapping method getNormalizedConcentrationCorrection.
/**
* TotalConservationCorrection is the term that takes local units (micro-molar) to either volume normalized micro-molar or area normalized molecules/sq-um
* @return cbit.vcell.parser.Expression
*/
public Expression getNormalizedConcentrationCorrection(SimulationContext simulationContext, UnitFactorProvider unitFactorProvider) throws ExpressionException, MappingException {
if (getGeometryClass() instanceof CompartmentSubVolume) {
if (simulationContext.getGeometryContext().isAllSizeSpecifiedPositive()) {
//
// everything mapped to micro-molar : just need size
//
Expression exp = new Expression(getSizeParameter(), simulationContext.getNameScope());
return exp;
} else {
throw new MappingException("\nIn non-spatial application '" + simulationContext.getName() + "', " + "size of structure '" + getStructure().getName() + "' must be assigned a " + "positive value if referenced in the model.\n\nPlease go to 'Structure Mapping' tab to check the size.");
}
} else if (getGeometryClass() instanceof SubVolume) {
//
// everything mapped to micro-molar : just need volume fraction
//
Expression exp = new Expression(getVolumePerUnitVolumeParameter(), simulationContext.getNameScope());
return exp;
} else if (getGeometryClass() instanceof SurfaceClass) {
//
// everything mapped to molecules/sq-um : need volume/area fraction and KMOLE
//
Expression unitFactor = unitFactorProvider.getUnitFactor(modelUnitSystem.getMembraneSubstanceUnit().divideBy(modelUnitSystem.getVolumeSubstanceUnit()));
Expression exp = Expression.mult(new Expression(getVolumePerUnitAreaParameter(), simulationContext.getNameScope()), unitFactor);
return exp;
} else {
throw new RuntimeException("structure " + getStructure().getName() + " not mapped");
}
}
Also used :
CompartmentSubVolume(cbit.vcell.geometry.CompartmentSubVolume)
Expression(cbit.vcell.parser.Expression)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
SubVolume(cbit.vcell.geometry.SubVolume)
CompartmentSubVolume(cbit.vcell.geometry.CompartmentSubVolume)
Example 18 with SurfaceClass
use of cbit.vcell.geometry.SurfaceClass in project vcell by virtualcell.
the class GeometryContext method refreshStructureMappings.
/**
* This method was created by a SmartGuide.
* @throws ExpressionBindingException
*/
public void refreshStructureMappings() throws MappingException, PropertyVetoException {
//
// step through all structure mappings
//
StructureMapping[] newStructureMappings = (StructureMapping[]) fieldStructureMappings.clone();
// fill in geometryClass for those MembraneMappings that were not mapped explicitly (and topology is available)
for (int j = 0; j < newStructureMappings.length; j++) {
StructureMapping structureMapping = newStructureMappings[j];
if (structureMapping instanceof MembraneMapping && structureMapping.getGeometryClass() == null) {
MembraneMapping membraneMapping = (MembraneMapping) structureMapping;
Feature insideFeature = getModel().getStructureTopology().getInsideFeature(membraneMapping.getMembrane());
Feature outsideFeature = getModel().getStructureTopology().getOutsideFeature(membraneMapping.getMembrane());
if (insideFeature != null && outsideFeature != null) {
FeatureMapping insideFeatureMapping = (FeatureMapping) getStructureMapping(insideFeature);
FeatureMapping outsideFeatureMapping = (FeatureMapping) getStructureMapping(outsideFeature);
if (insideFeatureMapping.getGeometryClass() == outsideFeatureMapping.getGeometryClass()) {
membraneMapping.setGeometryClass(insideFeatureMapping.getGeometryClass());
} else if (insideFeatureMapping.getGeometryClass() instanceof SubVolume && outsideFeatureMapping.getGeometryClass() instanceof SubVolume) {
SubVolume insideSubVolume = (SubVolume) insideFeatureMapping.getGeometryClass();
SubVolume outsideSubVolume = (SubVolume) outsideFeatureMapping.getGeometryClass();
SurfaceClass surfaceClass = getGeometry().getGeometrySurfaceDescription().getSurfaceClass(insideSubVolume, outsideSubVolume);
if (surfaceClass != null) {
membraneMapping.setGeometryClass(surfaceClass);
}
}
}
}
}
for (int j = 0; j < newStructureMappings.length; j++) {
StructureMapping structureMapping = newStructureMappings[j];
Structure mappedStructure = structureMapping.getStructure();
// SubVolume mappedSubvolume = structureMapping.getSubVolume();
Structure newStructure = null;
GeometryClass newGeometryClass = null;
boolean structureFound = false;
boolean geometryClassFound = false;
//
// match up with structures defined within model
//
Structure[] modelStructures = getModel().getStructures();
for (int i = 0; i < modelStructures.length; i++) {
Structure modelStructure = modelStructures[i];
if (modelStructure.compareEqual(mappedStructure)) {
structureFound = true;
newStructure = modelStructure;
break;
}
}
//
// match up with geometryClasses defined within geometry
//
GeometryClass[] geometryClasses = getGeometry().getGeometryClasses();
for (int i = 0; i < geometryClasses.length; i++) {
if (geometryClasses[i].compareEqual(structureMapping.getGeometryClass())) {
geometryClassFound = true;
newGeometryClass = geometryClasses[i];
break;
}
}
//
if (!(structureFound && geometryClassFound)) {
newStructureMappings = (StructureMapping[]) BeanUtils.removeElement(newStructureMappings, structureMapping);
j--;
// //
// // delete accompanied membrane mapping if exists
// //
// for (int i = 0; i < newStructureMappings.length; i++){
// if (newStructureMappings[i] instanceof MembraneMapping){
// MembraneMapping membraneMapping = (MembraneMapping)newStructureMappings[i];
// if (membraneMapping.getMembrane()==null ||
// membraneMapping.getMembrane().getInsideFeature() == structureMapping.getStructure() ||
// membraneMapping.getMembrane().getOutsideFeature() == structureMapping.getStructure()){
// newStructureMappings = (StructureMapping[])BeanUtils.removeElement(newStructureMappings,membraneMapping);
// break;
// }
// }
// }
} else {
// update references to Structure and SubVolume to correspond to those of Model and Geometry
structureMapping.setGeometryClass(newGeometryClass);
}
if (structureFound) {
structureMapping.setStructure(newStructure);
}
}
//
// add default mappings for any new structures
//
Structure[] structures = getModel().getStructures();
for (int i = 0; i < structures.length; i++) {
Structure structure = structures[i];
StructureMapping sm = null;
for (int j = 0; j < newStructureMappings.length; j++) {
if (newStructureMappings[j].getStructure().compareEqual(structure)) {
sm = newStructureMappings[j];
}
}
if (sm == null) {
if (structure instanceof Feature) {
FeatureMapping fm = new FeatureMapping((Feature) structure, fieldSimulationContext, getModel().getUnitSystem());
fm.setSimulationContext(this.fieldSimulationContext);
newStructureMappings = (StructureMapping[]) BeanUtils.addElement(newStructureMappings, fm);
if (getGeometry().getDimension() == 0) {
fm.setGeometryClass((CompartmentSubVolume) getGeometry().getGeometrySpec().getSubVolumes()[0]);
}
} else if (structure instanceof Membrane) {
MembraneMapping mm = new MembraneMapping((Membrane) structure, fieldSimulationContext, getModel().getUnitSystem());
mm.setSimulationContext(fieldSimulationContext);
newStructureMappings = (StructureMapping[]) BeanUtils.addElement(newStructureMappings, mm);
if (getGeometry().getDimension() == 0) {
mm.setGeometryClass((CompartmentSubVolume) getGeometry().getGeometrySpec().getSubVolumes()[0]);
}
} else {
throw new MappingException("unsupported Structure Mapping for structure " + structure.getClass().toString());
}
}
}
if (newStructureMappings != fieldStructureMappings) {
try {
setStructureMappings(newStructureMappings);
} catch (Exception e) {
e.printStackTrace(System.out);
throw new MappingException(e.getMessage());
}
}
fixMembraneMappings();
}
Also used :
GeometryClass(cbit.vcell.geometry.GeometryClass)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
Feature(cbit.vcell.model.Feature)
PropertyVetoException(java.beans.PropertyVetoException)
ExpressionException(cbit.vcell.parser.ExpressionException)
ExpressionBindingException(cbit.vcell.parser.ExpressionBindingException)
CompartmentSubVolume(cbit.vcell.geometry.CompartmentSubVolume)
SubVolume(cbit.vcell.geometry.SubVolume)
CompartmentSubVolume(cbit.vcell.geometry.CompartmentSubVolume)
Membrane(cbit.vcell.model.Membrane)
Structure(cbit.vcell.model.Structure)
Example 19 with SurfaceClass
use of cbit.vcell.geometry.SurfaceClass in project vcell by virtualcell.
the class GeometryContext method fixMembraneMappings.
private void fixMembraneMappings() throws PropertyVetoException {
StructureTopology structTopology = getModel().getStructureTopology();
for (int j = 0; j < fieldStructureMappings.length; j++) {
if (fieldStructureMappings[j] instanceof MembraneMapping) {
MembraneMapping membraneMapping = (MembraneMapping) fieldStructureMappings[j];
Membrane membrane = membraneMapping.getMembrane();
Feature insideFeature = structTopology.getInsideFeature(membrane);
Feature outsideFeature = structTopology.getOutsideFeature(membrane);
//
if (insideFeature != null && outsideFeature != null) {
FeatureMapping insideFM = (FeatureMapping) getStructureMapping(insideFeature);
FeatureMapping outsideFM = (FeatureMapping) getStructureMapping(outsideFeature);
GeometryClass insideGeometryClass = insideFM.getGeometryClass();
GeometryClass outsideGeometryClass = outsideFM.getGeometryClass();
//
if (insideFM != null && insideGeometryClass != null && outsideFM != null && outsideGeometryClass != null) {
// inside/outside both mapped to same domain ... membrane must be there too.
if (insideGeometryClass == outsideGeometryClass) {
membraneMapping.setGeometryClass(insideGeometryClass);
// inside/outside mapped to different subvolumes (try to map membrane to adjacent surfaceClass)
} else if (insideGeometryClass instanceof SubVolume && outsideGeometryClass instanceof SubVolume) {
GeometryClass[] geometryClasses = getGeometry().getGeometryClasses();
boolean bFound = false;
for (int i = 0; i < geometryClasses.length; i++) {
if (geometryClasses[i] instanceof SurfaceClass) {
SurfaceClass surfaceClass = (SurfaceClass) geometryClasses[i];
if (surfaceClass.isAdjacentTo((SubVolume) insideGeometryClass) && surfaceClass.isAdjacentTo((SubVolume) outsideGeometryClass)) {
membraneMapping.setGeometryClass(surfaceClass);
bFound = true;
}
}
}
if (!bFound) {
membraneMapping.setGeometryClass(null);
}
// inside/outside mapped to different membranes (membrane cannot be mapped ... must be cleared).
} else if (insideGeometryClass instanceof SurfaceClass && outsideGeometryClass instanceof SurfaceClass) {
membraneMapping.setGeometryClass(null);
// inside mapped to surface and outside mapped to subvolume (if adjacent, map membrane to surface ... else clear).
} else if (insideGeometryClass instanceof SurfaceClass && outsideGeometryClass instanceof SubVolume) {
SurfaceClass surface = (SurfaceClass) insideGeometryClass;
SubVolume subVolume = (SubVolume) outsideGeometryClass;
if (surface.isAdjacentTo(subVolume)) {
membraneMapping.setGeometryClass(surface);
} else {
membraneMapping.setGeometryClass(null);
}
// inside mapped to subvolume and outside mapped to surface (if adjacent, map membrane to surface ... else clear).
} else if (insideGeometryClass instanceof SubVolume && outsideGeometryClass instanceof SurfaceClass) {
SurfaceClass surface = (SurfaceClass) outsideGeometryClass;
SubVolume subVolume = (SubVolume) insideGeometryClass;
if (surface.isAdjacentTo(subVolume)) {
membraneMapping.setGeometryClass(surface);
} else {
membraneMapping.setGeometryClass(null);
}
}
}
}
}
}
}
Also used :
GeometryClass(cbit.vcell.geometry.GeometryClass)
StructureTopology(cbit.vcell.model.Model.StructureTopology)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
SubVolume(cbit.vcell.geometry.SubVolume)
CompartmentSubVolume(cbit.vcell.geometry.CompartmentSubVolume)
Membrane(cbit.vcell.model.Membrane)
Feature(cbit.vcell.model.Feature)
Example 20 with SurfaceClass
use of cbit.vcell.geometry.SurfaceClass in project vcell by virtualcell.
the class ParticleMathMapping method refreshMathDescription.
/**
* This method was created in VisualAge.
*/
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
getSimulationContext().checkValidity();
if (getSimulationContext().getGeometry().getDimension() == 0) {
throw new MappingException("particle math mapping requires spatial geometry - dimension >= 1");
}
StructureMapping[] structureMappings = getSimulationContext().getGeometryContext().getStructureMappings();
for (int i = 0; i < structureMappings.length; i++) {
if (structureMappings[i] instanceof MembraneMapping) {
if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
throw new MappingException("electric potential not yet supported for particle models");
}
}
}
//
// fail if any events
//
BioEvent[] bioEvents = getSimulationContext().getBioEvents();
if (bioEvents != null && bioEvents.length > 0) {
throw new MappingException("events not yet supported for particle-based models");
}
//
// gather only those reactionSteps that are not "excluded"
//
ReactionSpec[] reactionSpecs = getSimulationContext().getReactionContext().getReactionSpecs();
Vector<ReactionStep> rsList = new Vector<ReactionStep>();
for (int i = 0; i < reactionSpecs.length; i++) {
if (reactionSpecs[i].isExcluded() == false) {
if (reactionSpecs[i].isFast()) {
throw new MappingException("fast reactions not supported for particle models");
}
rsList.add(reactionSpecs[i].getReactionStep());
}
}
ReactionStep[] reactionSteps = new ReactionStep[rsList.size()];
rsList.copyInto(reactionSteps);
//
for (int i = 0; i < reactionSteps.length; i++) {
Kinetics.UnresolvedParameter[] unresolvedParameters = reactionSteps[i].getKinetics().getUnresolvedParameters();
if (unresolvedParameters != null && unresolvedParameters.length > 0) {
StringBuffer buffer = new StringBuffer();
for (int j = 0; j < unresolvedParameters.length; j++) {
if (j > 0) {
buffer.append(", ");
}
buffer.append(unresolvedParameters[j].getName());
}
throw new MappingException(reactionSteps[i].getDisplayType() + " '" + reactionSteps[i].getName() + "' contains unresolved identifier(s): " + buffer);
}
}
//
// temporarily place all variables in a hashtable (before binding) and discarding duplicates (check for equality)
//
VariableHash varHash = new VariableHash();
// //
// // verify that all structures are mapped to geometry classes and all geometry classes are mapped to a structure
// //
// Structure structures[] = getSimulationContext().getGeometryContext().getModel().getStructures();
// for (int i = 0; i < structures.length; i++){
// StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(structures[i]);
// if (sm==null || (sm.getGeometryClass() == null)){
// throw new MappingException("model structure '"+structures[i].getName()+"' not mapped to a geometry subdomain");
// }
// if (sm.getUnitSizeParameter()!=null){
// Expression unitSizeExp = sm.getUnitSizeParameter().getExpression();
// if(unitSizeExp != null)
// {
// try {
// double unitSize = unitSizeExp.evaluateConstant();
// if (unitSize != 1.0){
// throw new MappingException("model structure '"+sm.getStructure().getName()+"' unit size = "+unitSize+" != 1.0 ... partial volume or surface mapping not yet supported for particles");
// }
// }catch (ExpressionException e){
// e.printStackTrace(System.out);
// throw new MappingException("couldn't evaluate unit size for model structure '"+sm.getStructure().getName()+"' : "+e.getMessage());
// }
// }
// }
// }
// {
// GeometryClass[] geometryClass = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
// for (int i = 0; i < geometryClass.length; i++){
// Structure[] mappedStructures = getSimulationContext().getGeometryContext().getStructuresFromGeometryClass(geometryClass[i]);
// if (mappedStructures==null || mappedStructures.length==0){
// throw new MappingException("geometryClass '"+geometryClass[i].getName()+"' not mapped from a model structure");
// }
// }
// }
// deals with model parameters
Model model = getSimulationContext().getModel();
ModelUnitSystem modelUnitSystem = model.getUnitSystem();
ModelParameter[] modelParameters = model.getModelParameters();
// populate in globalParameterVariants hashtable
for (int j = 0; j < modelParameters.length; j++) {
Expression modelParamExpr = modelParameters[j].getExpression();
GeometryClass geometryClass = getDefaultGeometryClass(modelParamExpr);
modelParamExpr = getIdentifierSubstitutions(modelParamExpr, modelParameters[j].getUnitDefinition(), geometryClass);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), modelParamExpr, geometryClass));
}
//
// create new MathDescription (based on simContext's previous MathDescription if possible)
//
MathDescription oldMathDesc = getSimulationContext().getMathDescription();
mathDesc = null;
if (oldMathDesc != null) {
if (oldMathDesc.getVersion() != null) {
mathDesc = new MathDescription(oldMathDesc.getVersion());
} else {
mathDesc = new MathDescription(oldMathDesc.getName());
}
} else {
mathDesc = new MathDescription(getSimulationContext().getName() + "_generated");
}
//
// volume particle variables
//
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = enum1.nextElement();
if (scm.getVariable() instanceof ParticleVariable) {
if (!(mathDesc.getVariable(scm.getVariable().getName()) instanceof ParticleVariable)) {
varHash.addVariable(scm.getVariable());
}
}
}
varHash.addVariable(new Constant(getMathSymbol(model.getPI_CONSTANT(), null), getIdentifierSubstitutions(model.getPI_CONSTANT().getExpression(), model.getPI_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT().getExpression(), model.getFARADAY_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT_NMOLE(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT_NMOLE().getExpression(), model.getFARADAY_CONSTANT_NMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getGAS_CONSTANT(), null), getIdentifierSubstitutions(model.getGAS_CONSTANT().getExpression(), model.getGAS_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getTEMPERATURE(), null), getIdentifierSubstitutions(new Expression(getSimulationContext().getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
//
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping membraneMapping = (MembraneMapping) structureMappings[j];
GeometryClass geometryClass = membraneMapping.getGeometryClass();
//
// don't calculate voltage, still may need it though
//
Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
Variable voltageFunction = newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), geometryClass), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), geometryClass), geometryClass);
varHash.addVariable(voltageFunction);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), membraneMapping.getGeometryClass()), getIdentifierSubstitutions(membraneMapping.getInitialVoltageParameter().getExpression(), membraneMapping.getInitialVoltageParameter().getUnitDefinition(), membraneMapping.getGeometryClass()), membraneMapping.getGeometryClass()));
}
}
//
for (int j = 0; j < reactionSteps.length; j++) {
ReactionStep rs = reactionSteps[j];
if (getSimulationContext().getReactionContext().getReactionSpec(rs).isExcluded()) {
continue;
}
Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
GeometryClass geometryClass = null;
if (rs.getStructure() != null) {
geometryClass = getSimulationContext().getGeometryContext().getStructureMapping(rs.getStructure()).getGeometryClass();
}
if (parameters != null) {
for (int i = 0; i < parameters.length; i++) {
// Reaction rate, currentDensity, LumpedCurrent and null parameters are not going to displayed in the particle math description.
if (((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) || (parameters[i].getRole() == Kinetics.ROLE_LumpedCurrent) || (parameters[i].getRole() == Kinetics.ROLE_ReactionRate)) || (parameters[i].getExpression() == null)) {
continue;
}
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], geometryClass), getIdentifierSubstitutions(parameters[i].getExpression(), parameters[i].getUnitDefinition(), geometryClass), geometryClass));
}
}
}
//
// initial constants (either function or constant)
//
SpeciesContextSpec[] speciesContextSpecs = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpecParameter initParm = null;
Expression initExpr = null;
if (getSimulationContext().isUsingConcentration()) {
initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
initExpr = new Expression(initParm.getExpression());
// if (speciesContextSpecs[i].getSpeciesContext().getStructure() instanceof Feature) {
// initExpr = Expression.div(initExpr, new Expression(model.getKMOLE, getNameScope())).flatten();
// }
} else {
initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
initExpr = new Expression(initParm.getExpression());
}
if (initExpr != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
String[] symbols = initExpr.getSymbols();
// Check if 'initExpr' has other speciesContexts in its expression, need to replace it with 'spContext_init'
for (int j = 0; symbols != null && j < symbols.length; j++) {
// if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
SpeciesContext spC = null;
SymbolTableEntry ste = initExpr.getSymbolBinding(symbols[j]);
if (ste instanceof SpeciesContextSpecProxyParameter) {
SpeciesContextSpecProxyParameter spspp = (SpeciesContextSpecProxyParameter) ste;
if (spspp.getTarget() instanceof SpeciesContext) {
spC = (SpeciesContext) spspp.getTarget();
SpeciesContextSpec spcspec = getSimulationContext().getReactionContext().getSpeciesContextSpec(spC);
SpeciesContextSpecParameter spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
// if initConc param expression is null, try initCount
if (spCInitParm.getExpression() == null) {
spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
}
// need to get init condn expression, but can't get it from getMathSymbol() (mapping between bio and math), hence get it as below.
Expression scsInitExpr = new Expression(spCInitParm, getNameScope());
// scsInitExpr.bindExpression(this);
initExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
}
}
}
// now create the appropriate function for the current speciesContextSpec.
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParm, sm.getGeometryClass()), getIdentifierSubstitutions(initExpr, initParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter diffParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_DiffusionRate);
if (diffParm != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(diffParm, sm.getGeometryClass()), getIdentifierSubstitutions(diffParm.getExpression(), diffParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter bc_xm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXm);
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
if (bc_xm != null && (bc_xm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xm.getExpression(), bc_xm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_xp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXp);
if (bc_xp != null && (bc_xp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xp.getExpression(), bc_xp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_ym = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYm);
if (bc_ym != null && (bc_ym.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_ym, sm.getGeometryClass()), getIdentifierSubstitutions(bc_ym.getExpression(), bc_ym.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_yp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYp);
if (bc_yp != null && (bc_yp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_yp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_yp.getExpression(), bc_yp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZm);
if (bc_zm != null && (bc_zm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zm.getExpression(), bc_zm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZp);
if (bc_zp != null && (bc_zp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zp.getExpression(), bc_zp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter advection_velX = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityX);
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
GeometryClass geometryClass = sm.getGeometryClass();
if (advection_velX != null && (advection_velX.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velX, geometryClass), getIdentifierSubstitutions(advection_velX.getExpression(), advection_velX.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velY = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityY);
if (advection_velY != null && (advection_velY.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velY, geometryClass), getIdentifierSubstitutions(advection_velY.getExpression(), advection_velY.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velZ = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityZ);
if (advection_velZ != null && (advection_velZ.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velZ, geometryClass), getIdentifierSubstitutions(advection_velZ.getExpression(), advection_velZ.getUnitDefinition(), geometryClass), geometryClass));
}
}
//
// constant species (either function or constant)
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() instanceof Constant) {
varHash.addVariable(scm.getVariable());
}
}
//
// conversion factors
//
varHash.addVariable(new Constant(getMathSymbol(model.getKMOLE(), null), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getN_PMOLE(), null), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getKMILLIVOLTS(), null), getIdentifierSubstitutions(model.getKMILLIVOLTS().getExpression(), model.getKMILLIVOLTS().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getK_GHK(), null), getIdentifierSubstitutions(model.getK_GHK().getExpression(), model.getK_GHK().getUnitDefinition(), null)));
//
for (int i = 0; i < structureMappings.length; i++) {
StructureMapping sm = structureMappings[i];
if (getSimulationContext().getGeometry().getDimension() == 0) {
StructureMappingParameter sizeParm = sm.getSizeParameter();
if (sizeParm != null && sizeParm.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(sizeParm.getExpression(), sizeParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
} else {
if (sm instanceof MembraneMapping) {
MembraneMapping mm = (MembraneMapping) sm;
StructureMappingParameter volFrac = mm.getVolumeFractionParameter();
if (volFrac != null && volFrac.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(volFrac, sm.getGeometryClass()), getIdentifierSubstitutions(volFrac.getExpression(), volFrac.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
StructureMappingParameter surfToVol = mm.getSurfaceToVolumeParameter();
if (surfToVol != null && surfToVol.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(surfToVol, sm.getGeometryClass()), getIdentifierSubstitutions(surfToVol.getExpression(), surfToVol.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
} else {
Parameter parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
//
// functions
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
Variable dependentVariable = newFunctionOrConstant(getMathSymbol(scm.getSpeciesContext(), sm.getGeometryClass()), getIdentifierSubstitutions(scm.getDependencyExpression(), scm.getSpeciesContext().getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass());
dependentVariable.setDomain(new Domain(sm.getGeometryClass()));
varHash.addVariable(dependentVariable);
}
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass()));
}
}
//
// set Variables to MathDescription all at once with the order resolved by "VariableHash"
//
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
//
if (getSimulationContext().getGeometryContext().getGeometry() != null) {
try {
mathDesc.setGeometry(getSimulationContext().getGeometryContext().getGeometry());
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException("failure setting geometry " + e.getMessage());
}
} else {
throw new MappingException("geometry must be defined");
}
//
// create subdomains (volume and surfaces)
//
GeometryClass[] geometryClasses = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
for (int k = 0; k < geometryClasses.length; k++) {
if (geometryClasses[k] instanceof SubVolume) {
SubVolume subVolume = (SubVolume) geometryClasses[k];
//
// get priority of subDomain
//
// now does not have to match spatial feature, *BUT* needs to be unique
int priority = k;
//
// create subDomain
//
CompartmentSubDomain subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
mathDesc.addSubDomain(subDomain);
//
// assign boundary condition types
//
StructureMapping[] mappedSMs = getSimulationContext().getGeometryContext().getStructureMappings(subVolume);
FeatureMapping mappedFM = null;
for (int i = 0; i < mappedSMs.length; i++) {
if (mappedSMs[i] instanceof FeatureMapping) {
if (mappedFM != null) {
lg.warn("WARNING:::: MathMapping.refreshMathDescription() ... assigning boundary condition types not unique");
}
mappedFM = (FeatureMapping) mappedSMs[i];
}
}
if (mappedFM != null) {
subDomain.setBoundaryConditionXm(mappedFM.getBoundaryConditionTypeXm());
subDomain.setBoundaryConditionXp(mappedFM.getBoundaryConditionTypeXp());
if (getSimulationContext().getGeometry().getDimension() > 1) {
subDomain.setBoundaryConditionYm(mappedFM.getBoundaryConditionTypeYm());
subDomain.setBoundaryConditionYp(mappedFM.getBoundaryConditionTypeYp());
}
if (getSimulationContext().getGeometry().getDimension() > 2) {
subDomain.setBoundaryConditionZm(mappedFM.getBoundaryConditionTypeZm());
subDomain.setBoundaryConditionZp(mappedFM.getBoundaryConditionTypeZp());
}
}
} else if (geometryClasses[k] instanceof SurfaceClass) {
SurfaceClass surfaceClass = (SurfaceClass) geometryClasses[k];
// determine membrane inside and outside subvolume
// this preserves backward compatibility so that membrane subdomain
// inside and outside correspond to structure hierarchy when present
Pair<SubVolume, SubVolume> ret = DiffEquMathMapping.computeBoundaryConditionSource(model, simContext, surfaceClass);
SubVolume innerSubVolume = ret.one;
SubVolume outerSubVolume = ret.two;
//
// create subDomain
//
CompartmentSubDomain outerCompartment = mathDesc.getCompartmentSubDomain(outerSubVolume.getName());
CompartmentSubDomain innerCompartment = mathDesc.getCompartmentSubDomain(innerSubVolume.getName());
MembraneSubDomain memSubDomain = new MembraneSubDomain(innerCompartment, outerCompartment, surfaceClass.getName());
mathDesc.addSubDomain(memSubDomain);
}
}
//
// create Particle Contexts for all Particle Variables
//
Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
Expression unitFactor = getUnitFactor(modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getVolumeSubstanceUnit()));
while (enumSCM.hasMoreElements()) {
SpeciesContextMapping scm = enumSCM.nextElement();
SpeciesContext sc = scm.getSpeciesContext();
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure());
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
if (scm.getVariable() instanceof ParticleVariable && scm.getDependencyExpression() == null) {
ParticleVariable particleVariable = (ParticleVariable) scm.getVariable();
//
// initial distribution of particles
//
ArrayList<ParticleInitialCondition> particleInitialConditions = new ArrayList<ParticleInitialCondition>();
ParticleInitialCondition pic = null;
if (getSimulationContext().isUsingConcentration()) {
Expression initialDistribution = scs.getInitialConcentrationParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialConcentrationParameter(), sm.getGeometryClass()));
if (particleVariable instanceof VolumeParticleVariable) {
initialDistribution = Expression.mult(initialDistribution, unitFactor);
}
pic = new ParticleInitialConditionConcentration(initialDistribution);
} else {
Expression initialCount = scs.getInitialCountParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialCountParameter(), sm.getGeometryClass()));
if (initialCount == null) {
throw new MappingException("initialCount not defined for speciesContext " + scs.getSpeciesContext().getName());
}
Expression locationX = new Expression("u");
Expression locationY = new Expression("u");
Expression locationZ = new Expression("u");
pic = new ParticleInitialConditionCount(initialCount, locationX, locationY, locationZ);
}
particleInitialConditions.add(pic);
//
// diffusion
//
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm.getGeometryClass()));
Expression driftXExp = null;
if (scs.getVelocityXParameter().getExpression() != null) {
driftXExp = new Expression(getMathSymbol(scs.getVelocityXParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velX_quantities = scs.getVelocityQuantities(QuantityComponent.X);
if (velX_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velX_quantities.length == 1 && numRegions == 1) {
driftXExp = new Expression(getMathSymbol(velX_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
Expression driftYExp = null;
if (scs.getVelocityYParameter().getExpression() != null) {
driftYExp = new Expression(getMathSymbol(scs.getVelocityYParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velY_quantities = scs.getVelocityQuantities(QuantityComponent.Y);
if (velY_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velY_quantities.length == 1 && numRegions == 1) {
driftYExp = new Expression(getMathSymbol(velY_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
Expression driftZExp = null;
if (scs.getVelocityZParameter().getExpression() != null) {
driftZExp = new Expression(getMathSymbol(scs.getVelocityZParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velZ_quantities = scs.getVelocityQuantities(QuantityComponent.Z);
if (velZ_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velZ_quantities.length == 1 && numRegions == 1) {
driftZExp = new Expression(getMathSymbol(velZ_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
ParticleProperties particleProperties = new ParticleProperties(particleVariable, diffusion, driftXExp, driftYExp, driftZExp, particleInitialConditions);
GeometryClass myGC = sm.getGeometryClass();
if (myGC == null) {
throw new MappingException("Application '" + getSimulationContext().getName() + "'\nGeometry->StructureMapping->(" + sm.getStructure().getTypeName() + ")'" + sm.getStructure().getName() + "' must be mapped to geometry domain.\n(see 'Problems' tab)");
}
SubDomain subDomain = mathDesc.getSubDomain(myGC.getName());
subDomain.addParticleProperties(particleProperties);
}
}
for (ReactionStep reactionStep : reactionSteps) {
Kinetics kinetics = reactionStep.getKinetics();
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(reactionStep.getStructure());
GeometryClass reactionStepGeometryClass = sm.getGeometryClass();
SubDomain subdomain = mathDesc.getSubDomain(reactionStepGeometryClass.getName());
KineticsParameter reactionRateParameter = null;
if (kinetics instanceof LumpedKinetics) {
reactionRateParameter = ((LumpedKinetics) kinetics).getLumpedReactionRateParameter();
} else {
reactionRateParameter = ((DistributedKinetics) kinetics).getReactionRateParameter();
}
// macroscopic_irreversible/Microscopic_irreversible for bimolecular membrane reactions. They will NOT go through MassAction solver.
if (kinetics.getKineticsDescription().equals(KineticsDescription.Macroscopic_irreversible) || kinetics.getKineticsDescription().equals(KineticsDescription.Microscopic_irreversible)) {
Expression radiusExp = getIdentifierSubstitutions(reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_Binding_Radius).getExpression(), modelUnitSystem.getBindingRadiusUnit(), reactionStepGeometryClass);
if (radiusExp != null) {
Expression expCopy = new Expression(radiusExp);
try {
MassActionSolver.substituteParameters(expCopy, true).evaluateConstant();
} catch (ExpressionException e) {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Problem in binding radius of " + reactionStep.getName() + ": '" + radiusExp.infix() + "', " + e.getMessage()));
}
} else {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Binding radius of " + reactionStep.getName() + " is null."));
}
List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
List<Action> forwardActions = new ArrayList<Action>();
for (ReactionParticipant rp : reactionStep.getReactionParticipants()) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
if (rp instanceof Reactant) {
reactantParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (radiusExp != null) {
forwardActions.add(Action.createDestroyAction(particle));
}
}
}
} else if (rp instanceof Product) {
productParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (radiusExp != null) {
forwardActions.add(Action.createCreateAction(particle));
}
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
JumpProcessRateDefinition bindingRadius = new InteractionRadius(radiusExp);
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (forwardActions.size() > 0) {
ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, bindingRadius, forwardActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(forwardProcess);
}
} else // other type of reactions
{
/* check the reaction rate law to see if we need to decompose a reaction(reversible) into two jump processes.
rate constants are important in calculating the probability rate.
for Mass Action, we use KForward and KReverse,
for General Kinetics we parse reaction rate J to see if it is in Mass Action form.
*/
Expression forwardRate = null;
Expression reverseRate = null;
// Using the MassActionFunction to write out the math description
MassActionSolver.MassActionFunction maFunc = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction) || kinetics.getKineticsDescription().equals(KineticsDescription.General) || kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
Expression rateExp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
Parameter forwardRateParameter = null;
Parameter reverseRateParameter = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
} else if (kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
}
maFunc = MassActionSolver.solveMassAction(forwardRateParameter, reverseRateParameter, rateExp, reactionStep);
if (maFunc.getForwardRate() == null && maFunc.getReverseRate() == null) {
throw new MappingException("Cannot generate stochastic math mapping for the reaction:" + reactionStep.getName() + "\nLooking for the rate function according to the form of k1*Reactant1^Stoir1*Reactant2^Stoir2...-k2*Product1^Stoip1*Product2^Stoip2.");
} else {
if (maFunc.getForwardRate() != null) {
forwardRate = maFunc.getForwardRate();
}
if (maFunc.getReverseRate() != null) {
reverseRate = maFunc.getReverseRate();
}
}
}
if (maFunc != null) {
// if the reaction has forward rate (Mass action,HMMs), or don't have either forward or reverse rate (some other rate laws--like general)
// we process it as forward reaction
List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
List<Action> forwardActions = new ArrayList<Action>();
List<Action> reverseActions = new ArrayList<Action>();
List<ReactionParticipant> reactants = maFunc.getReactants();
List<ReactionParticipant> products = maFunc.getProducts();
for (ReactionParticipant rp : reactants) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
reactantParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (forwardRate != null) {
forwardActions.add(Action.createDestroyAction(particle));
}
if (reverseRate != null) {
reverseActions.add(Action.createCreateAction(particle));
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
for (ReactionParticipant rp : products) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
productParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (forwardRate != null) {
forwardActions.add(Action.createCreateAction(particle));
}
if (reverseRate != null) {
reverseActions.add(Action.createDestroyAction(particle));
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
//
// There are two unit conversions required:
//
// 1) convert entire reaction rate from vcell reaction units to Smoldyn units (molecules/lengthunit^dim/timeunit)
// (where dim is 2 for membrane reactions and 3 for volume reactions)
//
// for forward rates:
// 2) convert each reactant from Smoldyn units (molecules/lengthunit^dim) to VCell units
// (where dim is 2 for membrane reactants and 3 for volume reactants)
//
// or
//
// for reverse rates:
// 2) convert each product from Smoldyn units (molecules/lengthunit^dim) to VCell units
// (where dim is 2 for membrane products and 3 for volume products)
//
RationalNumber reactionLocationDim = new RationalNumber(reactionStep.getStructure().getDimension());
VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
VCUnitDefinition smoldynReactionSizeUnit = modelUnitSystem.getLengthUnit().raiseTo(reactionLocationDim);
VCUnitDefinition smoldynSubstanceUnit = modelUnitSystem.getStochasticSubstanceUnit();
VCUnitDefinition smoldynReactionRateUnit = smoldynSubstanceUnit.divideBy(smoldynReactionSizeUnit).divideBy(timeUnit);
VCUnitDefinition vcellReactionRateUnit = reactionRateParameter.getUnitDefinition();
VCUnitDefinition reactionUnitFactor = smoldynReactionRateUnit.divideBy(vcellReactionRateUnit);
if (forwardRate != null) {
VCUnitDefinition smoldynReactantsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
// start with factor to translate entire reaction rate.
VCUnitDefinition forwardUnitFactor = reactionUnitFactor;
//
for (ReactionParticipant reactant : maFunc.getReactants()) {
VCUnitDefinition vcellReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(reactant.getSpeciesContext()).isForceContinuous();
VCUnitDefinition smoldynReactantUnit = null;
if (bForceContinuous) {
// reactant is continuous (vcell units)
smoldynReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
} else {
// reactant is a particle (smoldyn units)
RationalNumber reactantLocationDim = new RationalNumber(reactant.getStructure().getDimension());
VCUnitDefinition smoldynReactantSize = modelUnitSystem.getLengthUnit().raiseTo(reactantLocationDim);
smoldynReactantUnit = smoldynSubstanceUnit.divideBy(smoldynReactantSize);
}
// keep track of units of all reactants
smoldynReactantsUnit = smoldynReactantsUnit.multiplyBy(smoldynReactantUnit);
RationalNumber reactantStoichiometry = new RationalNumber(reactant.getStoichiometry());
VCUnitDefinition reactantUnitFactor = (vcellReactantUnit.divideBy(smoldynReactantUnit)).raiseTo(reactantStoichiometry);
// accumulate unit factors for all reactants
forwardUnitFactor = forwardUnitFactor.multiplyBy(reactantUnitFactor);
}
forwardRate = Expression.mult(forwardRate, getUnitFactor(forwardUnitFactor));
VCUnitDefinition smoldynExpectedForwardRateUnit = smoldynReactionRateUnit.divideBy(smoldynReactantsUnit);
// get probability
Expression exp = getIdentifierSubstitutions(forwardRate, smoldynExpectedForwardRateUnit, reactionStepGeometryClass).flatten();
JumpProcessRateDefinition partRateDef = new MacroscopicRateConstant(exp);
// create particle jump process
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (forwardActions.size() > 0) {
ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, partRateDef, forwardActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(forwardProcess);
}
}
// end of forward rate not null
if (reverseRate != null) {
VCUnitDefinition smoldynProductsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
// start with factor to translate entire reaction rate.
VCUnitDefinition reverseUnitFactor = reactionUnitFactor;
//
for (ReactionParticipant product : maFunc.getProducts()) {
VCUnitDefinition vcellProductUnit = product.getSpeciesContext().getUnitDefinition();
boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(product.getSpeciesContext()).isForceContinuous();
VCUnitDefinition smoldynProductUnit = null;
if (bForceContinuous) {
smoldynProductUnit = product.getSpeciesContext().getUnitDefinition();
} else {
RationalNumber productLocationDim = new RationalNumber(product.getStructure().getDimension());
VCUnitDefinition smoldynProductSize = modelUnitSystem.getLengthUnit().raiseTo(productLocationDim);
smoldynProductUnit = smoldynSubstanceUnit.divideBy(smoldynProductSize);
}
// keep track of units of all products
smoldynProductsUnit = smoldynProductsUnit.multiplyBy(smoldynProductUnit);
RationalNumber productStoichiometry = new RationalNumber(product.getStoichiometry());
VCUnitDefinition productUnitFactor = (vcellProductUnit.divideBy(smoldynProductUnit)).raiseTo(productStoichiometry);
// accumulate unit factors for all products
reverseUnitFactor = reverseUnitFactor.multiplyBy(productUnitFactor);
}
reverseRate = Expression.mult(reverseRate, getUnitFactor(reverseUnitFactor));
VCUnitDefinition smoldynExpectedReverseRateUnit = smoldynReactionRateUnit.divideBy(smoldynProductsUnit);
// get probability
Expression exp = getIdentifierSubstitutions(reverseRate, smoldynExpectedReverseRateUnit, reactionStepGeometryClass).flatten();
JumpProcessRateDefinition partProbRate = new MacroscopicRateConstant(exp);
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName() + "_reverse");
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (reverseActions.size() > 0) {
ParticleJumpProcess reverseProcess = new ParticleJumpProcess(jpName, productParticles, partProbRate, reverseActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(reverseProcess);
}
}
// end of reverse rate not null
}
// end of maFunc not null
}
// end of reaction step for loop
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
if (mathDesc.getVariable(variable.getName()) == null) {
mathDesc.addVariable(variable);
}
}
}
if (!mathDesc.isValid()) {
lg.warn(mathDesc.getVCML_database());
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
if (lg.isDebugEnabled()) {
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
System.out.println(mathDesc.getVCML());
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
}
Also used :
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
MathDescription(cbit.vcell.math.MathDescription)
ArrayList(java.util.ArrayList)
Product(cbit.vcell.model.Product)
SpeciesContext(cbit.vcell.model.SpeciesContext)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Reactant(cbit.vcell.model.Reactant)
ExpressionException(cbit.vcell.parser.ExpressionException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SubVolume(cbit.vcell.geometry.SubVolume)
Vector(java.util.Vector)
SpatialQuantity(cbit.vcell.mapping.spatial.SpatialObject.SpatialQuantity)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition)
InteractionRadius(cbit.vcell.math.InteractionRadius)
ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
ParticleInitialCondition(cbit.vcell.math.ParticleProperties.ParticleInitialCondition)
ReactionStep(cbit.vcell.model.ReactionStep)
ParticleProperties(cbit.vcell.math.ParticleProperties)
Kinetics(cbit.vcell.model.Kinetics)
DistributedKinetics(cbit.vcell.model.DistributedKinetics)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
Domain(cbit.vcell.math.Variable.Domain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
ReactionParticipant(cbit.vcell.model.ReactionParticipant)
GeometryClass(cbit.vcell.geometry.GeometryClass)
Action(cbit.vcell.math.Action)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
Variable(cbit.vcell.math.Variable)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
VariableHash(cbit.vcell.math.VariableHash)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
Constant(cbit.vcell.math.Constant)
SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
RationalNumber(ucar.units_vcell.RationalNumber)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
Pair(org.vcell.util.Pair)
ParticleInitialConditionConcentration(cbit.vcell.math.ParticleProperties.ParticleInitialConditionConcentration)
ProcessSymmetryFactor(cbit.vcell.math.ParticleJumpProcess.ProcessSymmetryFactor)
Expression(cbit.vcell.parser.Expression)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MathException(cbit.vcell.math.MathException)
Model(cbit.vcell.model.Model)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Parameter(cbit.vcell.model.Parameter)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
MassActionSolver(cbit.vcell.model.MassActionSolver)
ParticleInitialConditionCount(cbit.vcell.math.ParticleProperties.ParticleInitialConditionCount)
Aggregations
SurfaceClass (cbit.vcell.geometry.SurfaceClass)48
SubVolume (cbit.vcell.geometry.SubVolume)42
GeometryClass (cbit.vcell.geometry.GeometryClass)19
Expression (cbit.vcell.parser.Expression)19
Geometry (cbit.vcell.geometry.Geometry)15
Feature (cbit.vcell.model.Feature)13
MathDescription (cbit.vcell.math.MathDescription)12
Model (cbit.vcell.model.Model)12
SpeciesContext (cbit.vcell.model.SpeciesContext)12
VCImage (cbit.image.VCImage)11
ImageSubVolume (cbit.vcell.geometry.ImageSubVolume)11
CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)11
MembraneSubDomain (cbit.vcell.math.MembraneSubDomain)11
Membrane (cbit.vcell.model.Membrane)11
Structure (cbit.vcell.model.Structure)11
FeatureMapping (cbit.vcell.mapping.FeatureMapping)10
ArrayList (java.util.ArrayList)10
ImageException (cbit.image.ImageException)9
AnalyticSubVolume (cbit.vcell.geometry.AnalyticSubVolume)8
CompartmentSubVolume (cbit.vcell.geometry.CompartmentSubVolume)8
