Example 21 with GeometrySpec
use of cbit.vcell.geometry.GeometrySpec in project vcell by virtualcell.
the class SBMLExporter method addGeometry.
private void addGeometry() throws SbmlException {
SpatialModelPlugin mplugin = (SpatialModelPlugin) sbmlModel.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
// Creates a geometry object via SpatialModelPlugin object.
org.sbml.jsbml.ext.spatial.Geometry sbmlGeometry = mplugin.createGeometry();
sbmlGeometry.setCoordinateSystem(GeometryKind.cartesian);
sbmlGeometry.setSpatialId("vcell");
Geometry vcGeometry = getSelectedSimContext().getGeometry();
Model vcModel = getSelectedSimContext().getModel();
//
// list of CoordinateComponents : 1 if geometry is 1-d, 2 if geometry is 2-d, 3 if geometry is 3-d
//
int dimension = vcGeometry.getDimension();
Extent vcExtent = vcGeometry.getExtent();
Origin vcOrigin = vcGeometry.getOrigin();
// add x coordinate component
CoordinateComponent xComp = sbmlGeometry.createCoordinateComponent();
xComp.setSpatialId(vcModel.getX().getName());
xComp.setType(CoordinateKind.cartesianX);
final UnitDefinition sbmlUnitDef_length = getOrCreateSBMLUnit(vcModel.getUnitSystem().getLengthUnit());
xComp.setUnits(sbmlUnitDef_length);
Boundary minX = new Boundary();
xComp.setBoundaryMinimum(minX);
minX.setSpatialId("Xmin");
minX.setValue(vcOrigin.getX());
Boundary maxX = new Boundary();
xComp.setBoundaryMaximum(maxX);
maxX.setSpatialId("Xmax");
maxX.setValue(vcOrigin.getX() + (vcExtent.getX()));
org.sbml.jsbml.Parameter pX = sbmlModel.createParameter();
pX.setId(vcModel.getX().getName());
pX.setValue(0.0);
pX.setConstant(false);
pX.setUnits(sbmlUnitDef_length);
SpatialParameterPlugin spPluginPx = (SpatialParameterPlugin) pX.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
SpatialSymbolReference spSymRefPx = new SpatialSymbolReference();
spPluginPx.setParamType(spSymRefPx);
spSymRefPx.setSpatialRef(xComp.getSpatialId());
// add y coordinate component
if (dimension == 2 || dimension == 3) {
CoordinateComponent yComp = sbmlGeometry.createCoordinateComponent();
yComp.setSpatialId(vcModel.getY().getName());
yComp.setType(CoordinateKind.cartesianY);
yComp.setUnits(sbmlUnitDef_length);
Boundary minY = new Boundary();
yComp.setBoundaryMinimum(minY);
minY.setSpatialId("Ymin");
minY.setValue(vcOrigin.getY());
Boundary maxY = new Boundary();
yComp.setBoundaryMaximum(maxY);
maxY.setSpatialId("Ymax");
maxY.setValue(vcOrigin.getY() + (vcExtent.getY()));
org.sbml.jsbml.Parameter pY = sbmlModel.createParameter();
pY.setId(vcModel.getY().getName());
pY.setValue(0.0);
pY.setConstant(false);
pY.setUnits(sbmlUnitDef_length);
SpatialParameterPlugin spPluginPy = (SpatialParameterPlugin) pY.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
SpatialSymbolReference spSymRefPy = new SpatialSymbolReference();
spPluginPy.setParamType(spSymRefPy);
spSymRefPy.setSpatialRef(yComp.getSpatialId());
}
// add z coordinate component
if (dimension == 3) {
CoordinateComponent zComp = sbmlGeometry.createCoordinateComponent();
zComp.setSpatialId(vcModel.getZ().getName());
zComp.setType(CoordinateKind.cartesianZ);
zComp.setUnits(sbmlUnitDef_length);
Boundary minZ = new Boundary();
zComp.setBoundaryMinimum(minZ);
minZ.setSpatialId("Zmin");
minZ.setValue(vcOrigin.getZ());
Boundary maxZ = new Boundary();
zComp.setBoundaryMaximum(maxZ);
maxZ.setSpatialId("Zmax");
maxZ.setValue(vcOrigin.getZ() + (vcExtent.getZ()));
org.sbml.jsbml.Parameter pZ = sbmlModel.createParameter();
pZ.setId(vcModel.getZ().getName());
pZ.setValue(0.0);
pZ.setConstant(false);
pZ.setUnits(sbmlUnitDef_length);
SpatialParameterPlugin spPluginPz = (SpatialParameterPlugin) pZ.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
SpatialSymbolReference spSymRefPz = new SpatialSymbolReference();
spPluginPz.setParamType(spSymRefPz);
spSymRefPz.setSpatialRef(zComp.getSpatialId());
}
//
// list of compartmentMappings : VC structureMappings
//
GeometryContext vcGeoContext = getSelectedSimContext().getGeometryContext();
StructureMapping[] vcStrucMappings = vcGeoContext.getStructureMappings();
for (int i = 0; i < vcStrucMappings.length; i++) {
StructureMapping vcStructMapping = vcStrucMappings[i];
String structName = vcStructMapping.getStructure().getName();
Compartment comp = sbmlModel.getCompartment(TokenMangler.mangleToSName(structName));
SpatialCompartmentPlugin cplugin = (SpatialCompartmentPlugin) comp.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
GeometryClass gc = vcStructMapping.getGeometryClass();
if (!goodPointer(gc, GeometryClass.class, structName)) {
continue;
}
CompartmentMapping compMapping = new CompartmentMapping();
cplugin.setCompartmentMapping(compMapping);
String geomClassName = gc.getName();
String id = TokenMangler.mangleToSName(geomClassName + structName);
compMapping.setSpatialId(id);
compMapping.setDomainType(TokenMangler.mangleToSName(DOMAIN_TYPE_PREFIX + geomClassName));
try {
StructureMappingParameter usp = vcStructMapping.getUnitSizeParameter();
Expression e = usp.getExpression();
if (goodPointer(e, Expression.class, id)) {
compMapping.setUnitSize(e.evaluateConstant());
}
} catch (ExpressionException e) {
e.printStackTrace(System.out);
throw new RuntimeException("Unable to create compartment mapping for structureMapping '" + compMapping.getId() + "' : " + e.getMessage());
}
}
//
// list of domain types : subvolumes and surface classes from VC
//
boolean bAnyAnalyticSubvolumes = false;
boolean bAnyImageSubvolumes = false;
boolean bAnyCSGSubvolumes = false;
GeometryClass[] vcGeomClasses = vcGeometry.getGeometryClasses();
int numSubVols = 0;
for (int i = 0; i < vcGeomClasses.length; i++) {
DomainType domainType = sbmlGeometry.createDomainType();
domainType.setSpatialId(DOMAIN_TYPE_PREFIX + vcGeomClasses[i].getName());
if (vcGeomClasses[i] instanceof SubVolume) {
if (((SubVolume) vcGeomClasses[i]) instanceof AnalyticSubVolume) {
bAnyAnalyticSubvolumes = true;
} else if (((SubVolume) vcGeomClasses[i]) instanceof ImageSubVolume) {
bAnyImageSubvolumes = true;
} else if (((SubVolume) vcGeomClasses[i]) instanceof CSGObject) {
bAnyCSGSubvolumes = true;
}
domainType.setSpatialDimensions(3);
numSubVols++;
} else if (vcGeomClasses[i] instanceof SurfaceClass) {
domainType.setSpatialDimensions(2);
}
}
//
// list of domains, adjacent domains : from VC geometricRegions
//
GeometrySurfaceDescription vcGSD = vcGeometry.getGeometrySurfaceDescription();
if (vcGSD.getRegionImage() == null) {
try {
vcGSD.updateAll();
} catch (Exception e) {
e.printStackTrace(System.out);
throw new RuntimeException("Unable to generate region images for geometry");
}
}
GeometricRegion[] vcGeometricRegions = vcGSD.getGeometricRegions();
ISize sampleSize = vcGSD.getVolumeSampleSize();
int numX = sampleSize.getX();
int numY = sampleSize.getY();
int numZ = sampleSize.getZ();
double ox = vcOrigin.getX();
double oy = vcOrigin.getY();
double oz = vcOrigin.getZ();
RegionInfo[] regionInfos = vcGSD.getRegionImage().getRegionInfos();
for (int i = 0; i < vcGeometricRegions.length; i++) {
// domains
Domain domain = sbmlGeometry.createDomain();
domain.setSpatialId(vcGeometricRegions[i].getName());
if (vcGeometricRegions[i] instanceof VolumeGeometricRegion) {
domain.setDomainType(DOMAIN_TYPE_PREFIX + ((VolumeGeometricRegion) vcGeometricRegions[i]).getSubVolume().getName());
//
// get a list of interior points ... should probably use the distance map to find a point
// furthest inside (or several points associated with the morphological skeleton).
//
InteriorPoint interiorPt = domain.createInteriorPoint();
int regionID = ((VolumeGeometricRegion) vcGeometricRegions[i]).getRegionID();
boolean bFound = false;
int regInfoIndx = 0;
for (int j = 0; j < regionInfos.length; j++) {
regInfoIndx = j;
if (regionInfos[j].getRegionIndex() == regionID) {
int volIndx = 0;
for (int z = 0; z < numZ && !bFound; z++) {
for (int y = 0; y < numY && !bFound; y++) {
for (int x = 0; x < numX && !bFound; x++) {
if (regionInfos[j].isIndexInRegion(volIndx)) {
bFound = true;
double unit_z = (numZ > 1) ? ((double) z) / (numZ - 1) : 0.5;
double coordZ = oz + vcExtent.getZ() * unit_z;
double unit_y = (numY > 1) ? ((double) y) / (numY - 1) : 0.5;
double coordY = oy + vcExtent.getY() * unit_y;
double unit_x = (numX > 1) ? ((double) x) / (numX - 1) : 0.5;
double coordX = ox + vcExtent.getX() * unit_x;
interiorPt.setCoord1(coordX);
interiorPt.setCoord2(coordY);
interiorPt.setCoord3(coordZ);
}
volIndx++;
}
// end - for x
}
// end - for y
}
// end - for z
}
// end if
}
// end for regionInfos
if (!bFound) {
throw new RuntimeException("Unable to find interior point for region '" + regionInfos[regInfoIndx].toString());
}
} else if (vcGeometricRegions[i] instanceof SurfaceGeometricRegion) {
SurfaceGeometricRegion vcSurfaceGeomReg = (SurfaceGeometricRegion) vcGeometricRegions[i];
GeometricRegion geomRegion0 = vcSurfaceGeomReg.getAdjacentGeometricRegions()[0];
GeometricRegion geomRegion1 = vcSurfaceGeomReg.getAdjacentGeometricRegions()[1];
SurfaceClass surfaceClass = vcGSD.getSurfaceClass(((VolumeGeometricRegion) geomRegion0).getSubVolume(), ((VolumeGeometricRegion) geomRegion1).getSubVolume());
domain.setDomainType(DOMAIN_TYPE_PREFIX + surfaceClass.getName());
// adjacent domains : 2 adjacent domain objects for each surfaceClass in VC.
// adjacent domain 1
GeometricRegion adjGeomRegion0 = vcSurfaceGeomReg.getAdjacentGeometricRegions()[0];
GeometricRegion adjGeomRegion1 = vcSurfaceGeomReg.getAdjacentGeometricRegions()[1];
AdjacentDomains adjDomain = new AdjacentDomains();
adjDomain.setSpatialId(TokenMangler.mangleToSName(vcSurfaceGeomReg.getName() + "_" + adjGeomRegion0.getName()));
adjDomain.setDomain1(vcSurfaceGeomReg.getName());
adjDomain.setDomain2(adjGeomRegion0.getName());
sbmlGeometry.addAdjacentDomain(adjDomain);
// adj domain 2
adjDomain = new AdjacentDomains();
adjDomain.setSpatialId(TokenMangler.mangleToSName(vcSurfaceGeomReg.getName() + "_" + adjGeomRegion1.getName()));
adjDomain.setDomain1(vcSurfaceGeomReg.getName());
adjDomain.setDomain2(adjGeomRegion1.getName());
sbmlGeometry.addAdjacentDomain(adjDomain);
}
}
//
if (bAnyAnalyticSubvolumes && !bAnyImageSubvolumes && !bAnyCSGSubvolumes) {
AnalyticGeometry sbmlAnalyticGeomDefinition = sbmlGeometry.createAnalyticGeometry();
sbmlAnalyticGeomDefinition.setSpatialId(TokenMangler.mangleToSName("Analytic_" + vcGeometry.getName()));
sbmlAnalyticGeomDefinition.setIsActive(true);
for (int i = 0; i < vcGeomClasses.length; i++) {
if (vcGeomClasses[i] instanceof AnalyticSubVolume) {
AnalyticVolume analyticVol = sbmlAnalyticGeomDefinition.createAnalyticVolume();
analyticVol.setSpatialId(vcGeomClasses[i].getName());
analyticVol.setDomainType(DOMAIN_TYPE_PREFIX + vcGeomClasses[i].getName());
analyticVol.setFunctionType(FunctionKind.layered);
analyticVol.setOrdinal(numSubVols - (i + 1));
Expression expr = ((AnalyticSubVolume) vcGeomClasses[i]).getExpression();
try {
String mathMLStr = ExpressionMathMLPrinter.getMathML(expr, true, MathType.BOOLEAN);
ASTNode mathMLNode = ASTNode.readMathMLFromString(mathMLStr);
analyticVol.setMath(mathMLNode);
} catch (Exception e) {
e.printStackTrace(System.out);
throw new RuntimeException("Error converting VC subvolume expression to mathML" + e.getMessage());
}
}
}
}
//
if (!bAnyAnalyticSubvolumes && !bAnyImageSubvolumes && bAnyCSGSubvolumes) {
CSGeometry sbmlCSGeomDefinition = new CSGeometry();
sbmlGeometry.addGeometryDefinition(sbmlCSGeomDefinition);
sbmlCSGeomDefinition.setSpatialId(TokenMangler.mangleToSName("CSG_" + vcGeometry.getName()));
for (int i = 0; i < vcGeomClasses.length; i++) {
if (vcGeomClasses[i] instanceof CSGObject) {
CSGObject vcellCSGObject = (CSGObject) vcGeomClasses[i];
org.sbml.jsbml.ext.spatial.CSGObject sbmlCSGObject = new org.sbml.jsbml.ext.spatial.CSGObject();
sbmlCSGeomDefinition.addCSGObject(sbmlCSGObject);
sbmlCSGObject.setSpatialId(vcellCSGObject.getName());
sbmlCSGObject.setDomainType(DOMAIN_TYPE_PREFIX + vcellCSGObject.getName());
// the ordinal should the the least for the default/background subVolume
sbmlCSGObject.setOrdinal(numSubVols - (i + 1));
org.sbml.jsbml.ext.spatial.CSGNode sbmlcsgNode = getSBMLCSGNode(vcellCSGObject.getRoot());
sbmlCSGObject.setCSGNode(sbmlcsgNode);
}
}
}
//
// add "Segmented" and "DistanceMap" SampledField Geometries
//
final boolean bVCGeometryIsImage = bAnyImageSubvolumes && !bAnyAnalyticSubvolumes && !bAnyCSGSubvolumes;
// 55if (bAnyAnalyticSubvolumes || bAnyImageSubvolumes || bAnyCSGSubvolumes){
if (bVCGeometryIsImage) {
//
// add "Segmented" SampledFieldGeometry
//
SampledFieldGeometry segmentedImageSampledFieldGeometry = sbmlGeometry.createSampledFieldGeometry();
segmentedImageSampledFieldGeometry.setSpatialId(TokenMangler.mangleToSName("SegmentedImage_" + vcGeometry.getName()));
segmentedImageSampledFieldGeometry.setIsActive(true);
// 55boolean bVCGeometryIsImage = bAnyImageSubvolumes && !bAnyAnalyticSubvolumes && !bAnyCSGSubvolumes;
Geometry vcImageGeometry = null;
{
if (bVCGeometryIsImage) {
// make a resampled image;
if (dimension == 3) {
try {
ISize imageSize = vcGeometry.getGeometrySpec().getDefaultSampledImageSize();
vcGeometry.precomputeAll(new GeometryThumbnailImageFactoryAWT());
vcImageGeometry = RayCaster.resampleGeometry(new GeometryThumbnailImageFactoryAWT(), vcGeometry, imageSize);
} catch (Throwable e) {
e.printStackTrace(System.out);
throw new RuntimeException("Unable to convert the original analytic or constructed solid geometry to image-based geometry : " + e.getMessage());
}
} else {
try {
vcGeometry.precomputeAll(new GeometryThumbnailImageFactoryAWT(), true, false);
GeometrySpec origGeometrySpec = vcGeometry.getGeometrySpec();
VCImage newVCImage = origGeometrySpec.getSampledImage().getCurrentValue();
//
// construct the new geometry with the sampled VCImage.
//
vcImageGeometry = new Geometry(vcGeometry.getName() + "_asImage", newVCImage);
vcImageGeometry.getGeometrySpec().setExtent(vcGeometry.getExtent());
vcImageGeometry.getGeometrySpec().setOrigin(vcGeometry.getOrigin());
vcImageGeometry.setDescription(vcGeometry.getDescription());
vcImageGeometry.getGeometrySurfaceDescription().setFilterCutoffFrequency(vcGeometry.getGeometrySurfaceDescription().getFilterCutoffFrequency());
vcImageGeometry.precomputeAll(new GeometryThumbnailImageFactoryAWT(), true, true);
} catch (Exception e) {
e.printStackTrace(System.out);
throw new RuntimeException("Unable to convert the original analytic or constructed solid geometry to image-based geometry : " + e.getMessage());
}
}
GeometryClass[] vcImageGeomClasses = vcImageGeometry.getGeometryClasses();
for (int j = 0; j < vcImageGeomClasses.length; j++) {
if (vcImageGeomClasses[j] instanceof ImageSubVolume) {
SampledVolume sampledVol = segmentedImageSampledFieldGeometry.createSampledVolume();
sampledVol.setSpatialId(vcGeomClasses[j].getName());
sampledVol.setDomainType(DOMAIN_TYPE_PREFIX + vcGeomClasses[j].getName());
sampledVol.setSampledValue(((ImageSubVolume) vcImageGeomClasses[j]).getPixelValue());
}
}
// add sampledField to sampledFieldGeometry
SampledField segmentedImageSampledField = sbmlGeometry.createSampledField();
VCImage vcImage = vcImageGeometry.getGeometrySpec().getImage();
segmentedImageSampledField.setSpatialId("SegmentedImageSampledField");
segmentedImageSampledField.setNumSamples1(vcImage.getNumX());
segmentedImageSampledField.setNumSamples2(vcImage.getNumY());
segmentedImageSampledField.setNumSamples3(vcImage.getNumZ());
segmentedImageSampledField.setInterpolationType(InterpolationKind.nearestneighbor);
segmentedImageSampledField.setCompression(CompressionKind.uncompressed);
segmentedImageSampledField.setDataType(DataKind.UINT8);
segmentedImageSampledFieldGeometry.setSampledField(segmentedImageSampledField.getId());
try {
byte[] vcImagePixelsBytes = vcImage.getPixels();
// imageData.setCompression("");
StringBuffer sb = new StringBuffer();
for (int i = 0; i < vcImagePixelsBytes.length; i++) {
int uint8_sample = ((int) vcImagePixelsBytes[i]) & 0xff;
sb.append(uint8_sample + " ");
}
segmentedImageSampledField.setSamplesLength(vcImage.getNumXYZ());
segmentedImageSampledField.setSamples(sb.toString().trim());
} catch (ImageException e) {
e.printStackTrace(System.out);
throw new RuntimeException("Unable to export image from VCell to SBML : " + e.getMessage());
}
}
}
/*
//
// add "DistanceMap" SampledFieldGeometry if there are exactly two subvolumes (else need more fields) and geometry is 3d.
//
if (numSubVols==2 && dimension == 3){
SignedDistanceMap[] distanceMaps = null;
try {
distanceMaps = DistanceMapGenerator.computeDistanceMaps(vcImageGeometry, vcImageGeometry.getGeometrySpec().getImage(), false, false);
} catch (ImageException e) {
e.printStackTrace(System.out);
System.err.println("Unable to export distance map sampled field from VCell to SBML : " + e.getMessage());
// throw new RuntimeException("Unable to export distance map sampled field from VCell to SBML : " + e.getMessage());
// don't want to throw an exception and stop export because distance map geometry couldn't be exported.
// just 'return' from method (since this is the last thing that is being done in this method).
return;
}
//
// the two distanceMaps should be redundant (one is negation of the other) ... so choose first one for field.
//
double[] signedDistances = distanceMaps[0].getSignedDistances();
SampledFieldGeometry distanceMapSampledFieldGeometry = sbmlGeometry.createSampledFieldGeometry();
distanceMapSampledFieldGeometry.setSpatialId(TokenMangler.mangleToSName("DistanceMap_"+vcGeometry.getName()));
SampledField distanceMapSampledField = distanceMapSampledFieldGeometry.createSampledField();
distanceMapSampledField.setSpatialId("DistanceMapSampledField");
distanceMapSampledField.setNumSamples1(distanceMaps[0].getSamplesX().length);
distanceMapSampledField.setNumSamples2(distanceMaps[0].getSamplesY().length);
distanceMapSampledField.setNumSamples3(distanceMaps[0].getSamplesZ().length);
distanceMapSampledField.setDataType("real");
System.err.println("do we need distanceMapSampleField.setDataType()?");
distanceMapSampledField.setInterpolationType("linear");
ImageData distanceMapImageData = distanceMapSampledField.createImageData();
distanceMapImageData.setDataType("int16");
System.err.println("should be:\n distanceMapImageData.setDataType(\"float32\")");
// distanceMapImageData.setCompression("");
double maxAbsValue = 0;
for (int i = 0; i < signedDistances.length; i++) {
maxAbsValue = Math.max(maxAbsValue,Math.abs(signedDistances[i]));
}
if (maxAbsValue==0.0){
throw new RuntimeException("computed distance map all zeros");
}
double scale = (Short.MAX_VALUE-1)/maxAbsValue;
int[] scaledIntegerDistanceMap = new int[signedDistances.length];
for (int i = 0; i < signedDistances.length; i++) {
scaledIntegerDistanceMap[i] = (int)(scale * signedDistances[i]);
}
distanceMapImageData.setSamples(scaledIntegerDistanceMap, signedDistances.length);
System.err.println("should be:\n distanceMapImageData.setSamples((float[])signedDistances,signedDistances.length)");
SampledVolume sampledVol = distanceMapSampledFieldGeometry.createSampledVolume();
sampledVol.setSpatialId(distanceMaps[0].getInsideSubvolumeName());
sampledVol.setDomainType(DOMAIN_TYPE_PREFIX+distanceMaps[0].getInsideSubvolumeName());
sampledVol.setSampledValue(255);
sampledVol = distanceMapSampledFieldGeometry.createSampledVolume();
sampledVol.setSpatialId(distanceMaps[1].getInsideSubvolumeName());
sampledVol.setDomainType(DOMAIN_TYPE_PREFIX+distanceMaps[1].getInsideSubvolumeName());
sampledVol.setSampledValue(1);
}
*/
}
//
// add "SurfaceMesh" ParametricGeometry
//
// if (bAnyAnalyticSubvolumes || bAnyImageSubvolumes || bAnyCSGSubvolumes){
// ParametricGeometry sbmlParametricGeomDefinition = sbmlGeometry.createParametricGeometry();
// sbmlParametricGeomDefinition.setSpatialId(TokenMangler.mangleToSName("SurfaceMesh_"+vcGeometry.getName()));
// xxxx
// }
}
Also used :
Origin(org.vcell.util.Origin)
CompartmentMapping(org.sbml.jsbml.ext.spatial.CompartmentMapping)
Compartment(org.sbml.jsbml.Compartment)
SpatialParameterPlugin(org.sbml.jsbml.ext.spatial.SpatialParameterPlugin)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
AnalyticGeometry(org.sbml.jsbml.ext.spatial.AnalyticGeometry)
ExpressionException(cbit.vcell.parser.ExpressionException)
Boundary(org.sbml.jsbml.ext.spatial.Boundary)
GeometrySpec(cbit.vcell.geometry.GeometrySpec)
DomainType(org.sbml.jsbml.ext.spatial.DomainType)
SampledVolume(org.sbml.jsbml.ext.spatial.SampledVolume)
SubVolume(cbit.vcell.geometry.SubVolume)
ImageSubVolume(cbit.vcell.geometry.ImageSubVolume)
AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume)
GeometryContext(cbit.vcell.mapping.GeometryContext)
SpatialCompartmentPlugin(org.sbml.jsbml.ext.spatial.SpatialCompartmentPlugin)
CoordinateComponent(org.sbml.jsbml.ext.spatial.CoordinateComponent)
SimulationContext(cbit.vcell.mapping.SimulationContext)
SpeciesContext(cbit.vcell.model.SpeciesContext)
GeometryContext(cbit.vcell.mapping.GeometryContext)
VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion)
SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion)
VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion)
GeometricRegion(cbit.vcell.geometry.surface.GeometricRegion)
SpatialSymbolReference(org.sbml.jsbml.ext.spatial.SpatialSymbolReference)
AnalyticVolume(org.sbml.jsbml.ext.spatial.AnalyticVolume)
InteriorPoint(org.sbml.jsbml.ext.spatial.InteriorPoint)
SampledField(org.sbml.jsbml.ext.spatial.SampledField)
Domain(org.sbml.jsbml.ext.spatial.Domain)
AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume)
GeometryClass(cbit.vcell.geometry.GeometryClass)
ImageException(cbit.image.ImageException)
GeometrySurfaceDescription(cbit.vcell.geometry.surface.GeometrySurfaceDescription)
Extent(org.vcell.util.Extent)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
ISize(org.vcell.util.ISize)
CSGeometry(org.sbml.jsbml.ext.spatial.CSGeometry)
RegionInfo(cbit.vcell.geometry.RegionImage.RegionInfo)
VCImage(cbit.image.VCImage)
StructureMapping(cbit.vcell.mapping.StructureMapping)
GeometryThumbnailImageFactoryAWT(cbit.vcell.geometry.GeometryThumbnailImageFactoryAWT)
ASTNode(org.sbml.jsbml.ASTNode)
CSGObject(cbit.vcell.geometry.CSGObject)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
UnitDefinition(org.sbml.jsbml.UnitDefinition)
ImageSubVolume(cbit.vcell.geometry.ImageSubVolume)
SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion)
SampledFieldGeometry(org.sbml.jsbml.ext.spatial.SampledFieldGeometry)
SpatialModelPlugin(org.sbml.jsbml.ext.spatial.SpatialModelPlugin)
InteriorPoint(org.sbml.jsbml.ext.spatial.InteriorPoint)
XMLStreamException(javax.xml.stream.XMLStreamException)
SbmlException(org.vcell.sbml.SbmlException)
ImageException(cbit.image.ImageException)
SBMLException(org.sbml.jsbml.SBMLException)
ExpressionException(cbit.vcell.parser.ExpressionException)
AdjacentDomains(org.sbml.jsbml.ext.spatial.AdjacentDomains)
Geometry(cbit.vcell.geometry.Geometry)
SampledFieldGeometry(org.sbml.jsbml.ext.spatial.SampledFieldGeometry)
AnalyticGeometry(org.sbml.jsbml.ext.spatial.AnalyticGeometry)
CSGeometry(org.sbml.jsbml.ext.spatial.CSGeometry)
Expression(cbit.vcell.parser.Expression)
Model(cbit.vcell.model.Model)
BioModel(cbit.vcell.biomodel.BioModel)
Example 22 with GeometrySpec
use of cbit.vcell.geometry.GeometrySpec in project vcell by virtualcell.
the class ComsolModelBuilder method getVCCModel.
public static VCCModel getVCCModel(SimulationJob vcellSimJob) throws ExpressionException {
MathDescription vcellMathDesc = vcellSimJob.getSimulation().getMathDescription();
Geometry vcellGeometry = vcellMathDesc.getGeometry();
GeometrySpec vcellGeometrySpec = vcellGeometry.getGeometrySpec();
int vcellDim = vcellGeometrySpec.getDimension();
VCCModel model = new VCCModel("Model", vcellDim);
model.modelpath = "D:\\Developer\\eclipse\\workspace_refactor\\comsol_java\\src";
model.comments = "Untitled\n\n";
VCCModelNode comp1 = new VCCModelNode("comp1");
model.modelnodes.add(comp1);
// if (vcellDim != 2){
// throw new RuntimeException("expecting 2D simulation");
// }
//
// assume initial geometry is circle centered at 0.5, 0.5 of radius 0.3
//
// String comsolOutsideDomainName = "dif1";
// String comsolInsideDomainName = "c1";
VCCGeomSequence geom1 = new VCCGeomSequence("geom1", vcellDim);
model.geometrysequences.add(geom1);
VCCMeshSequence mesh1 = new VCCMeshSequence("mesh1", geom1);
model.meshes.add(mesh1);
VCCStudy std1 = new VCCStudy("std1");
model.study = std1;
TimeBounds timeBounds = vcellSimJob.getSimulation().getSolverTaskDescription().getTimeBounds();
TimeStep timeStep = vcellSimJob.getSimulation().getSolverTaskDescription().getTimeStep();
String beginTime = Double.toString(timeBounds.getStartingTime());
String endTime = Double.toString(timeBounds.getEndingTime());
String step = Double.toString(timeStep.getDefaultTimeStep());
VCCStudyFeature time = new VCCTransientStudyFeature("time", beginTime, step, endTime);
std1.features.add(time);
if (vcellGeometrySpec.getImage() != null) {
throw new RuntimeException("image-based geometries not yet supported by VCell's COMSOL model builder");
}
if (vcellGeometrySpec.getNumSubVolumes() == 0) {
throw new RuntimeException("no subvolumes defined in geometry");
}
if (vcellGeometrySpec.getNumAnalyticOrCSGSubVolumes() != vcellGeometrySpec.getNumSubVolumes()) {
throw new RuntimeException("only analytic and CSG subvolumes currently supported by VCell's COMSOL model builder");
}
//
// add geometry for all subvolumes
//
HashMap<String, VCCGeomFeature> subvolumeNameFeatureMap = new HashMap<String, VCCGeomFeature>();
SubVolume[] subVolumes = vcellGeometrySpec.getSubVolumes();
for (int i = 0; i < subVolumes.length; i++) {
SubVolume subvolume = subVolumes[i];
if (subvolume instanceof CSGObject) {
CSGObject vcellCSGObject = (CSGObject) subvolume;
CSGNode vcellCSGNode = vcellCSGObject.getRoot();
ArrayList<VCCGeomFeature> geomFeatureList = new ArrayList<VCCGeomFeature>();
VCCGeomFeature feature = csgVisitor(vcellCSGNode, geomFeatureList, subvolume.getName());
geom1.geomfeatures.addAll(geomFeatureList);
if (i == 0) {
// first subvolume (on top in ordinals) doesn't need any differencing
subvolumeNameFeatureMap.put(subvolume.getName(), feature);
} else {
// have to subtract union of prior subvolumes
ArrayList<VCCGeomFeature> priorFeatures = new ArrayList<VCCGeomFeature>();
for (int j = 0; j < i; j++) {
CSGObject priorCSGObject = (CSGObject) subVolumes[j];
CSGNode priorCSGNode = priorCSGObject.getRoot();
geomFeatureList.clear();
VCCGeomFeature priorFeature = csgVisitor(priorCSGNode, geomFeatureList, subvolume.getName());
priorFeatures.add(priorFeature);
geom1.geomfeatures.addAll(geomFeatureList);
}
VCCDifference diff = new VCCDifference("diff" + subvolume.getName(), Keep.off);
diff.input.add(feature);
diff.input2.addAll(priorFeatures);
geom1.geomfeatures.add(diff);
subvolumeNameFeatureMap.put(subvolume.getName(), diff);
}
} else {
throw new RuntimeException("only CSG subvolumes currently supported by VCell's COMSOL model builder");
}
}
//
// add geometry for all surfaceClasses
//
HashMap<String, VCCGeomFeature> surfaceclassNameFeatureMap = new HashMap<String, VCCGeomFeature>();
SurfaceClass[] surfaceClasses = vcellGeometry.getGeometrySurfaceDescription().getSurfaceClasses();
for (int i = 0; i < surfaceClasses.length; i++) {
SurfaceClass surfaceClass = surfaceClasses[i];
Set<SubVolume> adjacentSubvolumes = surfaceClass.getAdjacentSubvolumes();
if (adjacentSubvolumes.size() != 2) {
throw new RuntimeException("expecting two adjacent subvolumes for surface " + surfaceClass.getName() + " in COMSOL model builder");
}
// find adjacent Geometry Features (for subvolumes)
Iterator<SubVolume> svIter = adjacentSubvolumes.iterator();
SubVolume subvolume0 = svIter.next();
SubVolume subvolume1 = svIter.next();
ArrayList<VCCGeomFeature> adjacentFeatures = new ArrayList<VCCGeomFeature>();
adjacentFeatures.add(subvolumeNameFeatureMap.get(subvolume0.getName()));
adjacentFeatures.add(subvolumeNameFeatureMap.get(subvolume1.getName()));
String name = "inter_" + subvolume0.getName() + "_" + subvolume1.getName();
// surfaces are dimension N-1
int entitydim = vcellDim - 1;
VCCIntersectionSelection intersect_subvolumes = new VCCIntersectionSelection(name, entitydim);
intersect_subvolumes.input.addAll(adjacentFeatures);
geom1.geomfeatures.add(intersect_subvolumes);
surfaceclassNameFeatureMap.put(surfaceClass.getName(), intersect_subvolumes);
}
SimulationSymbolTable symbolTable = new SimulationSymbolTable(vcellSimJob.getSimulation(), vcellSimJob.getJobIndex());
//
for (SubDomain subDomain : Collections.list(vcellMathDesc.getSubDomains())) {
for (Equation equ : subDomain.getEquationCollection()) {
if (equ instanceof PdeEquation || equ instanceof OdeEquation) {
VCCGeomFeature geomFeature = null;
final int dim;
if (subDomain instanceof CompartmentSubDomain) {
geomFeature = subvolumeNameFeatureMap.get(subDomain.getName());
dim = vcellDim;
} else if (subDomain instanceof MembraneSubDomain) {
geomFeature = surfaceclassNameFeatureMap.get(subDomain.getName());
dim = vcellDim - 1;
} else {
throw new RuntimeException("subdomains of type '" + subDomain.getClass().getSimpleName() + "' not yet supported in COMSOL model builder");
}
if (geomFeature == null) {
throw new RuntimeException("cannot find COMSOL geometry feature named " + subDomain.getName() + " in COMSOL model builder");
}
VCCConvectionDiffusionEquation cdeq = new VCCConvectionDiffusionEquation("cdeq_" + equ.getVariable().getName(), geom1, geomFeature, dim);
cdeq.fieldName = equ.getVariable().getName();
cdeq.initial = MathUtilities.substituteModelParameters(equ.getInitialExpression(), symbolTable).flatten().infix();
cdeq.sourceTerm_f = MathUtilities.substituteModelParameters(equ.getRateExpression(), symbolTable).flatten().infix();
if (equ instanceof PdeEquation) {
PdeEquation pde = (PdeEquation) equ;
cdeq.diffTerm_c = MathUtilities.substituteModelParameters(pde.getDiffusionExpression(), symbolTable).flatten().infix();
if (subDomain instanceof CompartmentSubDomain) {
CompartmentSubDomain compartmentSubdomain = (CompartmentSubDomain) subDomain;
ArrayList<String> be = new ArrayList<String>();
if (pde.getVelocityX() != null) {
be.add(MathUtilities.substituteModelParameters(pde.getVelocityX(), symbolTable).flatten().infix());
} else {
be.add("0");
}
if (vcellDim >= 2) {
if (pde.getVelocityY() != null) {
be.add(MathUtilities.substituteModelParameters(pde.getVelocityY(), symbolTable).flatten().infix());
} else {
be.add("0");
}
}
if (vcellDim == 3) {
if (pde.getVelocityY() != null) {
be.add(MathUtilities.substituteModelParameters(pde.getVelocityZ(), symbolTable).flatten().infix());
} else {
be.add("0");
}
}
cdeq.advection_be = be.toArray(new String[vcellDim]);
//
// look for membrane boundary conditions for this variable
//
MembraneSubDomain[] membraneSubdomains = vcellMathDesc.getMembraneSubDomains(compartmentSubdomain);
for (MembraneSubDomain membraneSubdomain : membraneSubdomains) {
JumpCondition jumpCondition = membraneSubdomain.getJumpCondition((VolVariable) pde.getVariable());
if (jumpCondition != null) {
Expression fluxExpr = null;
if (membraneSubdomain.getInsideCompartment() == compartmentSubdomain) {
fluxExpr = jumpCondition.getInFluxExpression();
} else if (membraneSubdomain.getOutsideCompartment() == compartmentSubdomain) {
fluxExpr = jumpCondition.getOutFluxExpression();
}
String name = equ.getVariable().getName() + "_flux_" + membraneSubdomain.getName();
VCCGeomFeature selection = surfaceclassNameFeatureMap.get(membraneSubdomain.getName());
VCCFluxBoundary fluxBoundary = new VCCFluxBoundary(name, selection, vcellDim - 1);
fluxBoundary.flux_g = MathUtilities.substituteModelParameters(fluxExpr, symbolTable).flatten().infix();
cdeq.features.add(fluxBoundary);
}
}
}
}
model.physics.add(cdeq);
}
}
}
//
return model;
}
Also used :
JumpCondition(cbit.vcell.math.JumpCondition)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
MathDescription(cbit.vcell.math.MathDescription)
HashMap(java.util.HashMap)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
VCCConvectionDiffusionEquation(org.vcell.solver.comsol.model.VCCConvectionDiffusionEquation)
CSGNode(cbit.vcell.geometry.CSGNode)
ArrayList(java.util.ArrayList)
VCCFluxBoundary(org.vcell.solver.comsol.model.VCCPhysicsFeature.VCCFluxBoundary)
GeometrySpec(cbit.vcell.geometry.GeometrySpec)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
PdeEquation(cbit.vcell.math.PdeEquation)
TimeBounds(cbit.vcell.solver.TimeBounds)
TimeStep(cbit.vcell.solver.TimeStep)
SubVolume(cbit.vcell.geometry.SubVolume)
VCCStudyFeature(org.vcell.solver.comsol.model.VCCStudyFeature)
CSGObject(cbit.vcell.geometry.CSGObject)
VCCGeomFeature(org.vcell.solver.comsol.model.VCCGeomFeature)
VCCDifference(org.vcell.solver.comsol.model.VCCGeomFeature.VCCDifference)
VCCGeomSequence(org.vcell.solver.comsol.model.VCCGeomSequence)
VCCIntersectionSelection(org.vcell.solver.comsol.model.VCCGeomFeature.VCCIntersectionSelection)
SimulationSymbolTable(cbit.vcell.solver.SimulationSymbolTable)
VCCConvectionDiffusionEquation(org.vcell.solver.comsol.model.VCCConvectionDiffusionEquation)
OdeEquation(cbit.vcell.math.OdeEquation)
PdeEquation(cbit.vcell.math.PdeEquation)
Equation(cbit.vcell.math.Equation)
VCCModelNode(org.vcell.solver.comsol.model.VCCModelNode)
VCCTransientStudyFeature(org.vcell.solver.comsol.model.VCCTransientStudyFeature)
VCCMeshSequence(org.vcell.solver.comsol.model.VCCMeshSequence)
Geometry(cbit.vcell.geometry.Geometry)
VCCStudy(org.vcell.solver.comsol.model.VCCStudy)
OdeEquation(cbit.vcell.math.OdeEquation)
Expression(cbit.vcell.parser.Expression)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
VCCModel(org.vcell.solver.comsol.model.VCCModel)
Example 23 with GeometrySpec
use of cbit.vcell.geometry.GeometrySpec in project vcell by virtualcell.
the class SmoldynFileWriter method writeWallSurfaces.
private void writeWallSurfaces() throws SolverException {
GeometrySurfaceDescription geometrySurfaceDescription = resampledGeometry.getGeometrySurfaceDescription();
GeometrySpec geometrySpec = resampledGeometry.getGeometrySpec();
SubVolume[] subVolumes = geometrySpec.getSubVolumes();
printWriter.println("# boundaries");
Origin origin = geometrySpec.getOrigin();
Extent extent = geometrySpec.getExtent();
Coordinate lowWall = new Coordinate(origin.getX(), origin.getY(), origin.getZ());
Coordinate highWall = new Coordinate(origin.getX() + extent.getX(), origin.getY() + extent.getY(), origin.getZ() + extent.getZ());
// potential artifact.
if (bHasNoSurface) {
SubDomain subDomain0 = mathDesc.getSubDomains().nextElement();
CompartmentSubDomain compartSubDomain0 = null;
compartSubDomain0 = (CompartmentSubDomain) subDomain0;
// x
if (compartSubDomain0.getBoundaryConditionXm().isPERIODIC()) {
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 0 " + lowWall.getX() + " " + highWall.getX() + " " + SmoldynVCellMapper.SmoldynKeyword.p);
} else {
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.low_wall + " 0 " + lowWall.getX() + " " + (compartSubDomain0.getBoundaryConditionXm().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.high_wall + " 0 " + highWall.getX() + " " + (compartSubDomain0.getBoundaryConditionXp().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
}
if (dimension > 1) {
// y
if (compartSubDomain0.getBoundaryConditionYm().isPERIODIC()) {
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 1 " + lowWall.getY() + " " + highWall.getY() + " " + SmoldynVCellMapper.SmoldynKeyword.p);
} else {
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.low_wall + " 1 " + lowWall.getY() + " " + (compartSubDomain0.getBoundaryConditionYm().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.high_wall + " 1 " + highWall.getY() + " " + (compartSubDomain0.getBoundaryConditionYp().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
}
if (dimension > 2) {
// z
if (compartSubDomain0.getBoundaryConditionZm().isPERIODIC()) {
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 2 " + lowWall.getZ() + " " + highWall.getZ() + " " + SmoldynVCellMapper.SmoldynKeyword.p);
} else {
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.low_wall + " 2 " + lowWall.getZ() + " " + (compartSubDomain0.getBoundaryConditionZm().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.high_wall + " 2 " + highWall.getZ() + " " + (compartSubDomain0.getBoundaryConditionZp().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
}
}
}
printWriter.println();
} else {
// x
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 0 " + lowWall.getX() + " " + highWall.getX());
if (dimension > 1) {
// y
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 1 " + lowWall.getY() + " " + highWall.getY());
if (dimension > 2) {
// z
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 2 " + lowWall.getZ() + " " + highWall.getZ());
}
}
printWriter.println();
// bounding walls as surfaces
// have to find boundary condition type
ISize sampleSize = simulation.getMeshSpecification().getSamplingSize();
int numX = sampleSize.getX();
int numY = dimension < 2 ? 1 : sampleSize.getY();
int numZ = dimension < 3 ? 1 : sampleSize.getZ();
if (dimension > 2) {
int[] k_wall = new int[] { 0, numZ - 1 };
for (int k = 0; k < k_wall.length; k++) {
for (int j = 0; j < numY; j++) {
for (int i = 0; i < numX; i++) {
int volIndex = k_wall[k] * numX * numY + j * numX + i;
for (SubVolume sv : subVolumes) {
// gather all the points in all the regions
GeometricRegion[] geometricRegions = geometrySurfaceDescription.getGeometricRegions(sv);
RegionInfo[] regionInfos = geometrySurfaceDescription.getRegionImage().getRegionInfos();
for (GeometricRegion gr : geometricRegions) {
VolumeGeometricRegion vgr = (VolumeGeometricRegion) gr;
for (RegionInfo ri : regionInfos) {
if (ri.getRegionIndex() == vgr.getRegionID() && ri.isIndexInRegion(volIndex)) {
boundaryZSubVolumes.add(sv);
}
}
}
}
}
}
}
}
if (dimension > 1) {
int[] j_wall = new int[] { 0, numY - 1 };
for (int k = 0; k < numZ; k++) {
for (int j = 0; j < j_wall.length; j++) {
for (int i = 0; i < numX; i++) {
int volIndex = k * numX * numY + j_wall[j] * numX + i;
for (SubVolume sv : subVolumes) {
// gather all the points in all the regions
GeometricRegion[] geometricRegions = geometrySurfaceDescription.getGeometricRegions(sv);
RegionInfo[] regionInfos = geometrySurfaceDescription.getRegionImage().getRegionInfos();
for (GeometricRegion gr : geometricRegions) {
VolumeGeometricRegion vgr = (VolumeGeometricRegion) gr;
for (RegionInfo ri : regionInfos) {
if (ri.getRegionIndex() == vgr.getRegionID() && ri.isIndexInRegion(volIndex)) {
boundaryYSubVolumes.add(sv);
}
}
}
}
}
}
}
}
int[] i_wall = new int[] { 0, numX - 1 };
for (int k = 0; k < numZ; k++) {
for (int j = 0; j < numY; j++) {
for (int i = 0; i < i_wall.length; i++) {
int volIndex = k * numX * numY + j * numX + i_wall[i];
for (SubVolume sv : subVolumes) {
// gather all the points in all the regions
GeometricRegion[] geometricRegions = geometrySurfaceDescription.getGeometricRegions(sv);
RegionInfo[] regionInfos = geometrySurfaceDescription.getRegionImage().getRegionInfos();
for (GeometricRegion gr : geometricRegions) {
VolumeGeometricRegion vgr = (VolumeGeometricRegion) gr;
for (RegionInfo ri : regionInfos) {
if (ri.getRegionIndex() == vgr.getRegionID() && ri.isIndexInRegion(volIndex)) {
boundaryXSubVolumes.add(sv);
}
}
}
}
}
}
}
Set<SubVolume> boundarySubVolumes = new HashSet<SubVolume>();
boundarySubVolumes.addAll(boundaryXSubVolumes);
boundarySubVolumes.addAll(boundaryYSubVolumes);
boundarySubVolumes.addAll(boundaryZSubVolumes);
BoundaryConditionType[] computedBct = new BoundaryConditionType[dimension * 2];
String[] smoldynBct = new String[dimension * 2];
String[] wallNames = new String[] { "Xm", "Xp", "Ym", "Yp", "Zm", "Zp" };
if (boundarySubVolumes.size() >= 1) {
for (SubVolume sv : boundarySubVolumes) {
CompartmentSubDomain csd = (CompartmentSubDomain) mathDesc.getSubDomain(sv.getName());
BoundaryConditionType[] bct = new BoundaryConditionType[] { csd.getBoundaryConditionXm(), csd.getBoundaryConditionXp(), csd.getBoundaryConditionYm(), csd.getBoundaryConditionYp(), csd.getBoundaryConditionZm(), csd.getBoundaryConditionZp() };
if (computedBct[0] == null) {
System.arraycopy(bct, 0, computedBct, 0, dimension * 2);
for (int i = 0; i < dimension * 2; i++) {
if (computedBct[i].isPERIODIC()) {
throw new SolverException("Models with both surfaces and periodic boundary conditions are not supported yet.");
}
smoldynBct[i] = computedBct[i].isDIRICHLET() ? SmoldynVCellMapper.SmoldynKeyword.absorb.name() : SmoldynVCellMapper.SmoldynKeyword.reflect.name();
}
} else {
for (int i = 0; i < dimension * 2; i++) {
if (!computedBct[i].compareEqual(bct[i])) {
throw new SolverException(wallNames[i] + " wall has different boundary conditions");
}
}
}
}
}
printWriter.println("# bounding wall surface");
// X walls
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.start_surface + " " + VCellSmoldynKeyword.bounding_wall_surface_X);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + "(" + SmoldynVCellMapper.SmoldynKeyword.up + ") " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.reflect);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.front + " " + smoldynBct[0]);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.back + " " + smoldynBct[1]);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.color + " " + SmoldynVCellMapper.SmoldynKeyword.both + " 1 1 1");
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_panels + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " 2");
// yz walls
switch(dimension) {
case 1:
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +0 " + lowWall.getX());
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -0 " + highWall.getX());
break;
case 2:
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +0 " + lowWall.getX() + " " + lowWall.getY() + " " + extent.getY());
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -0 " + highWall.getX() + " " + lowWall.getY() + " " + extent.getY());
break;
case 3:
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +0 " + lowWall.getX() + " " + lowWall.getY() + " " + lowWall.getZ() + " " + extent.getY() + " " + extent.getZ());
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -0 " + highWall.getX() + " " + lowWall.getY() + " " + lowWall.getZ() + " " + extent.getY() + " " + extent.getZ());
break;
}
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.end_surface);
printWriter.println();
if (dimension > 1) {
// Y walls
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.start_surface + " " + VCellSmoldynKeyword.bounding_wall_surface_Y);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + "(" + SmoldynVCellMapper.SmoldynKeyword.up + ") " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.reflect);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.front + " " + smoldynBct[2]);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.back + " " + smoldynBct[3]);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.color + " " + SmoldynVCellMapper.SmoldynKeyword.both + " 1 1 1");
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_panels + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " 2");
// xz walls
switch(dimension) {
case 2:
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +1 " + lowWall.getX() + " " + lowWall.getY() + " " + extent.getX());
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -1 " + lowWall.getX() + " " + highWall.getY() + " " + extent.getX());
break;
case 3:
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +1 " + lowWall.getX() + " " + lowWall.getY() + " " + lowWall.getZ() + " " + extent.getX() + " " + extent.getZ());
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -1 " + lowWall.getX() + " " + highWall.getY() + " " + lowWall.getZ() + " " + extent.getX() + " " + extent.getZ());
break;
}
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.end_surface);
printWriter.println();
if (dimension > 2) {
// Z walls
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.start_surface + " " + VCellSmoldynKeyword.bounding_wall_surface_Z);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + "(" + SmoldynVCellMapper.SmoldynKeyword.up + ") " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.reflect);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.front + " " + smoldynBct[4]);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.back + " " + smoldynBct[5]);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.color + " " + SmoldynVCellMapper.SmoldynKeyword.both + " 1 1 1");
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_panels + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " 2");
// xy walls
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +2 " + lowWall.getX() + " " + lowWall.getY() + " " + lowWall.getZ() + " " + extent.getX() + " " + extent.getY());
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -2 " + lowWall.getX() + " " + lowWall.getY() + " " + highWall.getZ() + " " + extent.getX() + " " + extent.getY());
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.end_surface);
printWriter.println();
}
}
}
}
Also used :
Origin(org.vcell.util.Origin)
GeometrySurfaceDescription(cbit.vcell.geometry.surface.GeometrySurfaceDescription)
Extent(org.vcell.util.Extent)
ISize(org.vcell.util.ISize)
RegionInfo(cbit.vcell.geometry.RegionImage.RegionInfo)
BoundaryConditionType(cbit.vcell.math.BoundaryConditionType)
VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion)
SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion)
VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion)
GeometricRegion(cbit.vcell.geometry.surface.GeometricRegion)
GeometrySpec(cbit.vcell.geometry.GeometrySpec)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
Coordinate(org.vcell.util.Coordinate)
SubVolume(cbit.vcell.geometry.SubVolume)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SolverException(cbit.vcell.solver.SolverException)
HashSet(java.util.HashSet)
Aggregations
GeometrySpec (cbit.vcell.geometry.GeometrySpec)23
SubVolume (cbit.vcell.geometry.SubVolume)12
VCImage (cbit.image.VCImage)7
SurfaceClass (cbit.vcell.geometry.SurfaceClass)6
GeometrySurfaceDescription (cbit.vcell.geometry.surface.GeometrySurfaceDescription)6
Expression (cbit.vcell.parser.Expression)6
ArrayList (java.util.ArrayList)6
AnalyticSubVolume (cbit.vcell.geometry.AnalyticSubVolume)5
Geometry (cbit.vcell.geometry.Geometry)5
GeometryThumbnailImageFactoryAWT (cbit.vcell.geometry.GeometryThumbnailImageFactoryAWT)5
GeometricRegion (cbit.vcell.geometry.surface.GeometricRegion)5
SurfaceGeometricRegion (cbit.vcell.geometry.surface.SurfaceGeometricRegion)5
VolumeGeometricRegion (cbit.vcell.geometry.surface.VolumeGeometricRegion)5
CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)5
MembraneSubDomain (cbit.vcell.math.MembraneSubDomain)5
Extent (org.vcell.util.Extent)5
Origin (org.vcell.util.Origin)5
CSGObject (cbit.vcell.geometry.CSGObject)4
ISize (org.vcell.util.ISize)4
VCImageUncompressed (cbit.image.VCImageUncompressed)3
