Example 1 with Node
use of cbit.vcell.geometry.surface.Node in project vcell by virtualcell.
the class FastMarchingMethodHA method main.
public static void main(String[] args) {
int numX = 2;
int numY = 2;
int numZ = 5;
int numItems = numX * numY * numZ;
double[] distanceMap = new double[numItems];
Arrays.fill(distanceMap, MAX_NUMBER);
// initial condition (exact distances from origin)
distanceMap[0] = 0;
distanceMap[1] = 1;
distanceMap[2] = 1;
// sqrt(2) - diagonal of a square
distanceMap[3] = 1.4142135;
distanceMap[4] = 1;
distanceMap[5] = 1.4142135;
distanceMap[6] = 1.4142135;
// sqrt(3) - diagonal of a cube
distanceMap[7] = 1.7320508;
// ----- next points distances for the above initial conditions
// FMM results vs. FMMHA results vs. exact computed results:
// 8: 2.000000 2.000000 2.000000
// 9: 2.350700 (0.12) 2.204818 (0.03) 2.236067 sqrt(5)
// 10: 2.350700 2.204818 2.236067
// 11: 2.642763 (0.20) 2.459807 (0.01) 2.449489 sqrt(6)
// 12: 3.000000 2.999999 3.000000
// 13: 3.303524 (0.14) 3.129413 (0.04) 3.162277 sqrt(10)
// 14: 3.303524 3.129413 3.162277 sqrt(10)
// 15: 3.569388 (0.25) 3.339079 (0.02) 3.316624 sqrt(2+9)
FastMarchingMethodHA fmm = new FastMarchingMethodHA(numX, numY, numZ, distanceMap, null);
System.out.println("Verifying indexes manipulation; index = x + y*numX + z*numX*numY");
int index = 0;
for (int z = 0; z < numZ; z++) {
for (int y = 0; y < numY; y++) {
for (int x = 0; x < numX; x++) {
index = x + y * numX + z * numX * numY;
System.out.println(x + ", " + y + ", " + z + ": " + index + " ");
System.out.println(fmm.getX(index) + ", " + fmm.getY(index) + ", " + fmm.getZ(index));
}
}
}
System.out.println("x, y, z: index");
fmm.march();
System.out.println("");
System.out.println(" -------------------- TEST 1 - accuracy test ----------------------- ");
long length = fmm.numX * fmm.numY * fmm.numZ;
for (index = 0; index < length; index++) {
System.out.println(index + ": " + distanceMap[index]);
}
System.out.println("");
System.out.println(" -------------------- TEST 2 - large array ----------------------- ");
numX = 100;
numY = 100;
numZ = 100;
numItems = numX * numY * numZ;
distanceMap = new double[numItems];
Arrays.fill(distanceMap, MAX_NUMBER);
// initial condition (a small cylinder) at 20, 20, 20
int x = 50;
int y = 50;
int z = 50;
int i = x + y * numX + z * numX * numY;
distanceMap[i] = 0;
distanceMap[i + 1] = 0;
distanceMap[i + 2] = 0;
// the tips
distanceMap[i - 1] = 1;
distanceMap[i + 3] = 1;
// closest neighbors on XZ plane
i = x + (y + 1) * numX + z * numX * numY;
distanceMap[i] = 0;
distanceMap[i + 1] = 1;
distanceMap[i + 2] = 1;
i = x + (y - 1) * numX + z * numX * numY;
distanceMap[i] = 0;
distanceMap[i + 1] = 1;
distanceMap[i + 2] = 1;
// closest neighbors on XY plane
i = x + y * numX + (z + 1) * numX * numY;
distanceMap[i] = 0;
distanceMap[i + 1] = 1;
distanceMap[i + 2] = 1;
i = x + y * numX + (z - 1) * numX * numY;
distanceMap[i] = 0;
distanceMap[i + 1] = 1;
distanceMap[i + 2] = 1;
fmm = new FastMarchingMethodHA(numX, numY, numZ, distanceMap, null);
long t1 = System.currentTimeMillis();
fmm.march();
long t2 = System.currentTimeMillis();
System.out.println("---------- " + (int) ((t2 - t1) / 1000) + " seconds ----------");
System.out.println("");
System.out.println(" -------------------- TEST 3 - compare accuracy ----------------------- ");
Vect3d p = null;
double d1 = 0;
numX = 100;
numY = 100;
numZ = 100;
numItems = numX * numY * numZ;
distanceMap = new double[numItems];
Arrays.fill(distanceMap, MAX_NUMBER);
Vect3d A = new Vect3d(0, 0, 0);
Vect3d B = new Vect3d(1, 0, 0);
Vect3d C = new Vect3d(0, 1, 0);
i = (int) (A.getX() + A.getY() * numX + A.getZ() * numX * numY);
distanceMap[i] = 0;
i = (int) (B.getX() + B.getY() * numX + B.getZ() * numX * numY);
distanceMap[i] = 0;
i = (int) (C.getX() + C.getY() * numX + C.getZ() * numX * numY);
distanceMap[i] = 0;
fmm = new FastMarchingMethodHA(numX, numY, numZ, distanceMap, null);
fmm.march();
final int numIterations = 500000;
double errorThreshold = 0.51;
int errorCount = 0;
Random rand = new Random();
try {
BufferedWriter out = new BufferedWriter(new FileWriter("c:\\TEMP\\FFM.3D"));
out.write("x y z value\n");
for (int j = 0; j < numIterations; j++) {
x = (int) (rand.nextDouble() * 100);
y = (int) (rand.nextDouble() * 100);
z = (int) (rand.nextDouble() * 100);
i = x + y * numX + z * numX * numY;
p = new Vect3d(x, y, z);
d1 = DistanceMapGenerator.distanceToTriangle3d(p, A, B, C);
double delta = Math.abs(d1 - distanceMap[i]);
String line = new String(x + " " + y + " " + z + " " + (int) (delta * 100) + "\n");
out.write(line);
if (delta > errorThreshold) {
System.out.println(x + ", " + y + ", " + z + " - delta: " + delta);
errorCount++;
}
}
out.close();
System.out.println("Delta above threshold " + errorThreshold + " occured in " + errorCount + " cases out of " + numIterations);
} catch (IOException e) {
e.printStackTrace();
}
System.out.println("A few random examples: ffm distance vs. exact distance");
for (int j = 0; j < 5; j++) {
x = (int) (1 + rand.nextDouble() * 98);
y = (int) (1 + rand.nextDouble() * 98);
z = (int) (1 + rand.nextDouble() * 98);
i = x + y * numX + z * numX * numY;
p = new Vect3d(x, y, z);
d1 = DistanceMapGenerator.distanceToTriangle3d(p, A, B, C);
System.out.println(x + ", " + y + ", " + z + " - ffm dist: " + distanceMap[i] + ", exact dist: " + d1);
}
System.out.println("");
System.out.println(" -------------------- TEST 4 - negative entries ----------------------- ");
p = null;
d1 = 0;
numX = 10;
numY = 10;
numZ = 10;
numItems = numX * numY * numZ;
distanceMap = new double[numItems];
Arrays.fill(distanceMap, MAX_NUMBER);
Node AA = new Node(0.501, 0.5, 0.5);
Node BB = new Node(0.5, 0.501, 0.5);
Node CC = new Node(0.5, 0.5, 0.501);
A = new Vect3d(AA);
B = new Vect3d(BB);
C = new Vect3d(CC);
distanceMap[0] = -0.86602540378443864676372317075294;
distanceMap[1] = 0.86602540378443864676372317075294;
distanceMap[10] = 0.86602540378443864676372317075294;
distanceMap[11] = 0.86602540378443864676372317075294;
distanceMap[100] = 0.86602540378443864676372317075294;
distanceMap[101] = 0.86602540378443864676372317075294;
distanceMap[110] = 0.86602540378443864676372317075294;
distanceMap[111] = 0.86602540378443864676372317075294;
fmm = new FastMarchingMethodHA(numX, numY, numZ, distanceMap, null);
fmm.march();
x = 5;
y = 5;
z = 5;
i = x + y * numX + z * numX * numY;
p = new Vect3d(x, y, z);
Node pp = new Node(x, y, z);
d1 = DistanceMapGenerator.distanceToTriangle3d(p, A, B, C);
double d2 = DistanceMapGenerator.distanceToTriangleExperimental(pp, AA, BB, CC, "c:\\TEMP\\triangle.3D");
System.out.println(x + ", " + y + ", " + z + " - ffm dist: " + distanceMap[i] + ", exact dist: " + d1 + ", exper dist: " + d2);
x = 2;
y = 1;
z = 0;
i = x + y * numX + z * numX * numY;
p = new Vect3d(x, y, z);
pp = new Node(x, y, z);
d1 = DistanceMapGenerator.distanceToTriangle3d(p, A, B, C);
d2 = DistanceMapGenerator.distanceToTriangleExperimental(pp, AA, BB, CC, "c:\\TEMP\\triangle.3D");
System.out.println(x + ", " + y + ", " + z + " - ffm dist: " + distanceMap[i] + ", exact dist: " + d1 + ", exper dist: " + d2);
System.out.println("");
System.out.println(" -------------------- TEST 5 ----------------------- ");
p = null;
d1 = 0;
numX = 10;
numY = 10;
numZ = 10;
numItems = numX * numY * numZ;
distanceMap = new double[numItems];
Arrays.fill(distanceMap, MAX_NUMBER);
AA = new Node(0.01, 0, 0);
BB = new Node(0, 0.01, 0);
CC = new Node(0, 0, 0.01);
A = new Vect3d(AA);
B = new Vect3d(BB);
C = new Vect3d(CC);
distanceMap[0] = 0;
distanceMap[1] = 1;
distanceMap[10] = 1;
distanceMap[11] = 1.4142135623730950488016887242097;
distanceMap[100] = 1;
distanceMap[101] = 1.4142135623730950488016887242097;
distanceMap[110] = 1.4142135623730950488016887242097;
distanceMap[111] = 1.7320508075688772935274463415059;
fmm = new FastMarchingMethodHA(numX, numY, numZ, distanceMap, null);
fmm.march();
x = 5;
y = 5;
z = 5;
i = x + y * numX + z * numX * numY;
p = new Vect3d(x, y, z);
pp = new Node(x, y, z);
d1 = DistanceMapGenerator.distanceToTriangle3d(p, A, B, C);
d2 = DistanceMapGenerator.distanceToTriangleExperimental(pp, AA, BB, CC, "c:\\TEMP\\triangle.3D");
System.out.println(x + ", " + y + ", " + z + " - ffm dist: " + distanceMap[i] + ", exact dist: " + d1 + ", exper dist: " + d2);
x = 2;
y = 1;
z = 0;
i = x + y * numX + z * numX * numY;
p = new Vect3d(x, y, z);
pp = new Node(x, y, z);
d1 = DistanceMapGenerator.distanceToTriangle3d(p, A, B, C);
d2 = DistanceMapGenerator.distanceToTriangleExperimental(pp, AA, BB, CC, "c:\\TEMP\\triangle.3D");
System.out.println(x + ", " + y + ", " + z + " - ffm dist: " + distanceMap[i] + ", exact dist: " + d1 + ", exper dist: " + d2);
System.out.println("");
System.out.println(" -------------------- TEST 6 - narrow vertical column ---------------------- ");
p = null;
d1 = 0;
numX = 2;
numY = 2;
numZ = 10;
numItems = numX * numY * numZ;
distanceMap = new double[numItems];
Arrays.fill(distanceMap, MAX_NUMBER);
distanceMap[0] = -0.3;
distanceMap[1] = -0.3;
distanceMap[2] = -0.3;
distanceMap[3] = -0.3;
distanceMap[4] = 0.7;
distanceMap[5] = 0.7;
distanceMap[6] = 0.7;
distanceMap[7] = 0.7;
fmm = new FastMarchingMethodHA(numX, numY, numZ, distanceMap, null);
fmm.march();
for (int j = 8; j < numItems; j++) {
System.out.println(j + ": " + distanceMap[j]);
}
System.out.println("");
System.out.println(" -------------------- TEST 7 - 2D test, numZ must be 1 ---------------------- ");
p = null;
d1 = 0;
numX = 500;
numY = 500;
numZ = 1;
numItems = numX * numY * numZ;
distanceMap = new double[numItems];
Arrays.fill(distanceMap, MAX_NUMBER);
x = 0;
y = 0;
double radius = 200;
for (int alpha = 0; alpha < 360; alpha++) {
double rad = Math.PI * (double) alpha / 180.0;
x = (int) (250.0 + radius * Math.cos(rad));
y = (int) (250.0 + radius * Math.sin(rad));
i = x + y * numX;
distanceMap[i] = 0;
}
try {
// save some points in a VisIt compatible format
BufferedWriter out = new BufferedWriter(new FileWriter("c:\\TEMP\\FFM2D1" + ".3D"));
out.write("x y z value\n");
for (int j = 0; j < distanceMap.length; j++) {
x = DistanceMapGenerator.getX(j, numX, numY);
y = DistanceMapGenerator.getY(j, numX, numY);
z = DistanceMapGenerator.getZ(j, numX, numY);
if (distanceMap[j] < MAX_NUMBER) {
out.write(x + " " + y + " " + z + " " + (int) (distanceMap[j] * 10) + "\n");
}
}
out.close();
} catch (IOException e) {
}
fmm = new FastMarchingMethodHA(numX, numY, numZ, distanceMap, null);
fmm.march();
try {
// save some points in a VisIt compatible format
BufferedWriter out = new BufferedWriter(new FileWriter("c:\\TEMP\\FFM2D2" + ".3D"));
out.write("x y z value\n");
for (int j = 0; j < distanceMap.length; j++) {
x = DistanceMapGenerator.getX(j, numX, numY);
y = DistanceMapGenerator.getY(j, numX, numY);
z = DistanceMapGenerator.getZ(j, numX, numY);
if (distanceMap[j] < MAX_NUMBER) {
out.write(x + " " + y + " " + z + " " + (int) (distanceMap[j] * 10) + "\n");
}
}
out.close();
} catch (IOException e) {
}
System.out.println("");
System.out.println(" -------------------- TEST 8 - anisotropic ---------------------- ");
numX = 2;
numY = 2;
numZ = 5;
numItems = numX * numY * numZ;
distanceMap = new double[numItems];
Arrays.fill(distanceMap, MAX_NUMBER);
// initial condition (exact distances from origin)
distanceMap[0] = 0;
distanceMap[1] = 1;
distanceMap[2] = 1;
// sqrt(2) - diagonal of a square
distanceMap[3] = 1.4142135623730950488016887242097;
distanceMap[4] = 2;
// sqrt(5)
distanceMap[5] = 2.2360679774997896964091736687313;
distanceMap[6] = 2.2360679774997896964091736687313;
// sqrt(6)
distanceMap[7] = 2.4494897427831780981972840747059;
// ----- next points distances for the above initial conditions
// FMM-plain FMM-HA exact computed:
// 8: 4.00000000 N/A 4.00000000
// 9: 4.19691079 4.12310562 // sqrt(17)
// 10:
// 11: 4.38072927 4.24264068
// 12: 6.00000000 6.00000000
// 13: 6.16836145 6.08276253 // sqrt(37)
// 14:
// 15: 6.32866014 6.16441400 // sqrt(16+2)
fmm = new FastMarchingMethodHA(numX, numY, numZ, 1, 1, 2, distanceMap, null);
System.out.println("Verifying indexes manipulation; index = x + y*numX + z*numX*numY");
fmm.march();
length = fmm.numX * fmm.numY * fmm.numZ;
for (index = 0; index < length; index++) {
System.out.println(index + ": " + distanceMap[index]);
}
System.out.println("");
System.out.println(" -------------------- TEST 9 - LONG test ---------------------- ");
numX = 300;
numY = 300;
numZ = 300;
numItems = numX * numY * numZ;
distanceMap = new double[numItems];
Arrays.fill(distanceMap, MAX_NUMBER);
distanceMap[0] = 0;
distanceMap[1] = 1;
distanceMap[300] = 1;
distanceMap[301] = 1.4142135623730950488016887242097;
distanceMap[90000] = 1;
distanceMap[90001] = 1.4142135623730950488016887242097;
distanceMap[90300] = 1.4142135623730950488016887242097;
distanceMap[90301] = 1.7320508075688772935274463415059;
fmm = new FastMarchingMethodHA(numX, numY, numZ, distanceMap, null);
fmm.march();
System.out.println("Done");
}
Also used :
Random(java.util.Random)
FileWriter(java.io.FileWriter)
Node(cbit.vcell.geometry.surface.Node)
IOException(java.io.IOException)
BufferedWriter(java.io.BufferedWriter)
Example 2 with Node
use of cbit.vcell.geometry.surface.Node in project vcell by virtualcell.
the class SmoldynFileWriter method writeSurfaces.
private void writeSurfaces() throws SolverException, ImageException, PropertyVetoException, GeometryException, ExpressionException {
GeometrySurfaceDescription geometrySurfaceDescription = resampledGeometry.getGeometrySurfaceDescription();
SurfaceClass[] surfaceClasses = geometrySurfaceDescription.getSurfaceClasses();
GeometrySpec geometrySpec = resampledGeometry.getGeometrySpec();
SubVolume[] surfaceGeometrySubVolumes = geometrySpec.getSubVolumes();
GeometricRegion[] AllGeometricRegions = resampledGeometry.getGeometrySurfaceDescription().getGeometricRegions();
ArrayList<SurfaceGeometricRegion> surfaceRegionList = new ArrayList<SurfaceGeometricRegion>();
ArrayList<VolumeGeometricRegion> volumeRegionList = new ArrayList<VolumeGeometricRegion>();
for (GeometricRegion geometricRegion : AllGeometricRegions) {
if (geometricRegion instanceof SurfaceGeometricRegion) {
surfaceRegionList.add((SurfaceGeometricRegion) geometricRegion);
} else if (geometricRegion instanceof VolumeGeometricRegion) {
volumeRegionList.add((VolumeGeometricRegion) geometricRegion);
} else {
throw new SolverException("unsupported geometric region type " + geometricRegion.getClass());
}
}
printWriter.println("# geometry");
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.dim + " " + dimension);
if (bHasNoSurface) {
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_compartment + " " + surfaceGeometrySubVolumes.length);
} else {
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_compartment + " " + (surfaceGeometrySubVolumes.length + 1));
// plus the surface which are bounding walls
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_surface + " " + (surfaceClasses.length + dimension));
}
printWriter.println();
// write boundaries and wall surfaces
writeWallSurfaces();
// for 3D ... smoldyn normal convension is triangle right-hand-rule normal points to the outside compartment subdomain.
if (!bHasNoSurface) {
membraneSubdomainTriangleMap = new HashMap<MembraneSubDomain, ArrayList<TrianglePanel>>();
// write surfaces
printWriter.println("# surfaces");
int triangleGlobalCount = 0;
int membraneIndex = -1;
SurfaceCollection surfaceCollection = geometrySurfaceDescription.getSurfaceCollection();
// pre-allocate collections used repeatedly in following loops; clear before reusing
HashMap<Node, Set<String>> nodeTriMap = new HashMap<>();
ArrayList<TrianglePanel> triList = new ArrayList<TrianglePanel>();
// use a sorted set to ensure neighbors written out is same order for reproducibility
SortedSet<String> neighborsForCurrentNode = new TreeSet<String>();
for (int sci = 0; sci < surfaceClasses.length; sci++) {
nodeTriMap.clear();
triList.clear();
int triLocalCount = 0;
SurfaceClass surfaceClass = surfaceClasses[sci];
GeometricRegion[] geometricRegions = geometrySurfaceDescription.getGeometricRegions(surfaceClass);
for (GeometricRegion gr : geometricRegions) {
SurfaceGeometricRegion sgr = (SurfaceGeometricRegion) gr;
VolumeGeometricRegion volRegion0 = (VolumeGeometricRegion) sgr.getAdjacentGeometricRegions()[0];
VolumeGeometricRegion volRegion1 = (VolumeGeometricRegion) sgr.getAdjacentGeometricRegions()[1];
SubVolume subVolume0 = volRegion0.getSubVolume();
SubVolume subVolume1 = volRegion1.getSubVolume();
CompartmentSubDomain compart0 = mathDesc.getCompartmentSubDomain(subVolume0.getName());
CompartmentSubDomain compart1 = mathDesc.getCompartmentSubDomain(subVolume1.getName());
MembraneSubDomain membraneSubDomain = mathDesc.getMembraneSubDomain(compart0, compart1);
if (membraneSubDomain == null) {
throw new SolverException(VCellErrorMessages.getSmoldynUnexpectedSurface(compart0, compart1));
}
int exteriorRegionID = volRegion0.getRegionID();
int interiorRegionID = volRegion1.getRegionID();
if (membraneSubDomain.getInsideCompartment() == compart0) {
exteriorRegionID = volRegion1.getRegionID();
interiorRegionID = volRegion0.getRegionID();
}
for (int j = 0; j < surfaceCollection.getSurfaceCount(); j++) {
Surface surface = surfaceCollection.getSurfaces(j);
if ((surface.getInteriorRegionIndex() == exteriorRegionID && surface.getExteriorRegionIndex() == interiorRegionID) || (surface.getInteriorRegionIndex() == interiorRegionID && surface.getExteriorRegionIndex() == exteriorRegionID)) {
// Polygon polygon = surface.getPolygons(k);
for (Polygon polygon : surface) {
if (polygonMembaneElementMap != null) {
membraneIndex = polygonMembaneElementMap.get(polygon).getMembraneIndex();
}
Node[] nodes = polygon.getNodes();
if (dimension == 2) {
// ignore z
Vect3d unitNormal = new Vect3d();
polygon.getUnitNormal(unitNormal);
unitNormal.set(unitNormal.getX(), unitNormal.getY(), 0);
int point0 = 0;
Vect3d v0 = new Vect3d(nodes[point0].getX(), nodes[point0].getY(), 0);
int point1 = 1;
Vect3d v1 = null;
for (point1 = 1; point1 < nodes.length; point1++) {
if (v0.getX() != nodes[point1].getX() || v0.getY() != nodes[point1].getY()) {
v1 = new Vect3d(nodes[point1].getX(), nodes[point1].getY(), 0);
break;
}
}
if (v1 == null) {
throw new RuntimeException("failed to generate surface");
}
Vect3d v01 = Vect3d.sub(v1, v0);
Vect3d unit01n = v01.cross(unitNormal);
unit01n.unit();
if (Math.abs(unit01n.getZ() - 1.0) < 1e-6) {
// v0 to v1 opposes vcell surface normal. it's already flipped.
Triangle triangle;
if (surface.getInteriorRegionIndex() == interiorRegionID) {
// we have to flipped it back
triangle = new Triangle(nodes[point1], nodes[point0], null);
} else {
triangle = new Triangle(nodes[point0], nodes[point1], null);
}
triList.add(new TrianglePanel(triLocalCount++, triangleGlobalCount++, membraneIndex, triangle));
} else if (Math.abs(unit01n.getZ() + 1.0) < 1e-6) {
// v0 to v1 is in direction of vcell surface normal.
Triangle triangle;
if (surface.getInteriorRegionIndex() == interiorRegionID) {
triangle = new Triangle(nodes[point0], nodes[point1], null);
} else {
triangle = new Triangle(nodes[point1], nodes[point0], null);
}
triList.add(new TrianglePanel(triLocalCount++, triangleGlobalCount++, membraneIndex, triangle));
} else {
throw new RuntimeException("failed to generate surface");
}
} else if (dimension == 3) {
Triangle triangle1;
Triangle triangle2;
if (surface.getInteriorRegionIndex() == interiorRegionID) {
// interior
triangle1 = new Triangle(nodes[0], nodes[1], nodes[2]);
triangle2 = new Triangle(nodes[0], nodes[2], nodes[3]);
} else {
triangle1 = new Triangle(nodes[2], nodes[1], nodes[0]);
triangle2 = new Triangle(nodes[3], nodes[2], nodes[0]);
}
triList.add(new TrianglePanel(triLocalCount++, triangleGlobalCount++, membraneIndex, triangle1));
triList.add(new TrianglePanel(triLocalCount++, triangleGlobalCount++, membraneIndex, triangle2));
}
}
}
}
}
// add triangles to node hash
for (TrianglePanel triPanel : triList) {
for (Node node : triPanel.triangle.getNodes()) {
if (node == null) {
continue;
}
Set<String> triNameSet = nodeTriMap.get(node);
if (triNameSet == null) {
triNameSet = new HashSet<String>();
nodeTriMap.put(node, triNameSet);
}
triNameSet.add(triPanel.name);
}
}
SubVolume[] adjacentSubvolums = surfaceClass.getAdjacentSubvolumes().toArray(new SubVolume[0]);
CompartmentSubDomain csd0 = simulation.getMathDescription().getCompartmentSubDomain(adjacentSubvolums[0].getName());
CompartmentSubDomain csd1 = simulation.getMathDescription().getCompartmentSubDomain(adjacentSubvolums[1].getName());
MembraneSubDomain membraneSubDomain = simulation.getMathDescription().getMembraneSubDomain(csd0, csd1);
membraneSubdomainTriangleMap.put(membraneSubDomain, triList);
final boolean initialMoleculesOnMembrane = (closestTriangles != null);
if (initialMoleculesOnMembrane) {
findClosestTriangles(membraneSubDomain, triList);
}
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.start_surface + " " + surfaceClass.getName());
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + "(" + SmoldynVCellMapper.SmoldynKeyword.all + ") " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.reflect);
// printWriter.println(SmoldynKeyword.action + " " + SmoldynKeyword.all + "(" + SmoldynKeyword.up + ") " + SmoldynKeyword.both + " " + SmoldynKeyword.reflect);
// get color after species
Color c = colors[sci + particleVariableList.size()];
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.color + " " + SmoldynVCellMapper.SmoldynKeyword.both + " " + c.getRed() / 255.0 + " " + c.getGreen() / 255.0 + " " + c.getBlue() / 255.0 + " 0.1");
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.front + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.back + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_panels + " " + SmoldynVCellMapper.SmoldynKeyword.tri + " " + triList.size());
for (TrianglePanel trianglePanel : triList) {
Triangle triangle = trianglePanel.triangle;
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.tri);
switch(dimension) {
case 1:
printWriter.print(" " + triangle.getNodes(0).getX());
break;
case 2:
printWriter.print(" " + triangle.getNodes(0).getX() + " " + triangle.getNodes(0).getY());
printWriter.print(" " + triangle.getNodes(1).getX() + " " + triangle.getNodes(1).getY());
break;
case 3:
for (Node node : triangle.getNodes()) {
printWriter.print(" " + node.getX() + " " + node.getY() + " " + node.getZ());
}
break;
}
printWriter.println(" " + trianglePanel.name);
}
for (TrianglePanel triPanel : triList) {
neighborsForCurrentNode.clear();
for (Node node : triPanel.triangle.getNodes()) {
if (node == null) {
continue;
}
neighborsForCurrentNode.addAll(nodeTriMap.get(node));
}
neighborsForCurrentNode.remove(triPanel.name);
// printWriter.print(SmoldynKeyword.neighbors + " " +triPanel.name);
// to allow smoldyn read line length as 256, chop the neighbors to multiple lines
int maxNeighborCount = 4;
//
int count = 0;
for (String neigh : neighborsForCurrentNode) {
if (count % maxNeighborCount == 0) {
printWriter.println();
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.neighbors + " " + triPanel.name);
}
printWriter.print(" " + neigh);
count++;
}
}
printWriter.println();
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.end_surface);
printWriter.println();
}
// write compartment
// printWriter.println("# bounding wall compartment");
// printWriter.println(SmoldynKeyword.start_compartment + " " + VCellSmoldynKeyword.bounding_wall_compartment);
// printWriter.println(SmoldynKeyword.surface + " " + VCellSmoldynKeyword.bounding_wall_surface_X);
// if (dimension > 1) {
// printWriter.println(SmoldynKeyword.surface + " " + VCellSmoldynKeyword.bounding_wall_surface_Y);
// if (dimension > 2) {
// printWriter.println(SmoldynKeyword.surface + " " + VCellSmoldynKeyword.bounding_wall_surface_Z);
// }
// }
// printWriter.println(SmoldynKeyword.end_compartment);
// printWriter.println();
}
}
Also used :
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
Set(java.util.Set)
TreeSet(java.util.TreeSet)
SortedSet(java.util.SortedSet)
DataSet(cbit.vcell.simdata.DataSet)
HashSet(java.util.HashSet)
GeometrySurfaceDescription(cbit.vcell.geometry.surface.GeometrySurfaceDescription)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
HashMap(java.util.HashMap)
Node(cbit.vcell.geometry.surface.Node)
ArrayList(java.util.ArrayList)
Triangle(cbit.vcell.geometry.surface.Triangle)
Surface(cbit.vcell.geometry.surface.Surface)
GeometrySpec(cbit.vcell.geometry.GeometrySpec)
SubVolume(cbit.vcell.geometry.SubVolume)
TreeSet(java.util.TreeSet)
Polygon(cbit.vcell.geometry.surface.Polygon)
SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion)
SurfaceCollection(cbit.vcell.geometry.surface.SurfaceCollection)
Color(java.awt.Color)
VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion)
SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion)
VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion)
GeometricRegion(cbit.vcell.geometry.surface.GeometricRegion)
Vect3d(cbit.vcell.render.Vect3d)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SolverException(cbit.vcell.solver.SolverException)
Example 3 with Node
use of cbit.vcell.geometry.surface.Node in project vcell by virtualcell.
the class SmoldynSurfaceTessellator method writeSurfaces.
protected void writeSurfaces() throws SolverException, ImageException, PropertyVetoException, GeometryException, ExpressionException {
GeometrySurfaceDescription geometrySurfaceDescription = resampledGeometry.getGeometrySurfaceDescription();
SurfaceClass[] surfaceClasses = geometrySurfaceDescription.getSurfaceClasses();
GeometrySpec geometrySpec = resampledGeometry.getGeometrySpec();
SubVolume[] surfaceGeometrySubVolumes = geometrySpec.getSubVolumes();
GeometricRegion[] AllGeometricRegions = resampledGeometry.getGeometrySurfaceDescription().getGeometricRegions();
ArrayList<SurfaceGeometricRegion> surfaceRegionList = new ArrayList<SurfaceGeometricRegion>();
ArrayList<VolumeGeometricRegion> volumeRegionList = new ArrayList<VolumeGeometricRegion>();
for (GeometricRegion geometricRegion : AllGeometricRegions) {
if (geometricRegion instanceof SurfaceGeometricRegion) {
surfaceRegionList.add((SurfaceGeometricRegion) geometricRegion);
} else if (geometricRegion instanceof VolumeGeometricRegion) {
volumeRegionList.add((VolumeGeometricRegion) geometricRegion);
} else {
throw new SolverException("unsupported geometric region type " + geometricRegion.getClass());
}
}
printWriter.println("# geometry");
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.dim + " " + dimension);
if (bHasNoSurface) {
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_compartment + " " + surfaceGeometrySubVolumes.length);
} else {
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_compartment + " " + (surfaceGeometrySubVolumes.length + 1));
// plus the surface which are bounding walls
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_surface + " " + (surfaceClasses.length + dimension));
}
printWriter.println();
// write boundaries and wall surfaces
writeWallSurfaces();
// for 3D ... smoldyn normal convension is triangle right-hand-rule normal points to the outside compartment subdomain.
if (!bHasNoSurface) {
membraneSubdomainTriangleMap = new HashMap<MembraneSubDomain, ArrayList<TrianglePanel>>();
// write surfaces
printWriter.println("# surfaces");
int triangleGlobalCount = 0;
int membraneIndex = -1;
SurfaceCollection surfaceCollection = geometrySurfaceDescription.getSurfaceCollection();
// pre-allocate collections used repeatedly in following loops; clear before reusing
HashMap<Node, Set<String>> nodeTriMap = new HashMap<>();
ArrayList<TrianglePanel> triList = new ArrayList<TrianglePanel>();
// use a sorted set to ensure neighbors written out is same order for reproducibility
SortedSet<String> neighborsForCurrentNode = new TreeSet<String>();
for (int sci = 0; sci < surfaceClasses.length; sci++) {
nodeTriMap.clear();
triList.clear();
int triLocalCount = 0;
SurfaceClass surfaceClass = surfaceClasses[sci];
GeometricRegion[] geometricRegions = geometrySurfaceDescription.getGeometricRegions(surfaceClass);
for (GeometricRegion gr : geometricRegions) {
SurfaceGeometricRegion sgr = (SurfaceGeometricRegion) gr;
VolumeGeometricRegion volRegion0 = (VolumeGeometricRegion) sgr.getAdjacentGeometricRegions()[0];
VolumeGeometricRegion volRegion1 = (VolumeGeometricRegion) sgr.getAdjacentGeometricRegions()[1];
SubVolume subVolume0 = volRegion0.getSubVolume();
SubVolume subVolume1 = volRegion1.getSubVolume();
CompartmentSubDomain compart0 = mathDesc.getCompartmentSubDomain(subVolume0.getName());
CompartmentSubDomain compart1 = mathDesc.getCompartmentSubDomain(subVolume1.getName());
MembraneSubDomain membraneSubDomain = mathDesc.getMembraneSubDomain(compart0, compart1);
if (membraneSubDomain == null) {
throw new SolverException(VCellErrorMessages.getSmoldynUnexpectedSurface(compart0, compart1));
}
int exteriorRegionID = volRegion0.getRegionID();
int interiorRegionID = volRegion1.getRegionID();
if (membraneSubDomain.getInsideCompartment() == compart0) {
exteriorRegionID = volRegion1.getRegionID();
interiorRegionID = volRegion0.getRegionID();
}
for (int j = 0; j < surfaceCollection.getSurfaceCount(); j++) {
Surface surface = surfaceCollection.getSurfaces(j);
if ((surface.getInteriorRegionIndex() == exteriorRegionID && surface.getExteriorRegionIndex() == interiorRegionID) || (surface.getInteriorRegionIndex() == interiorRegionID && surface.getExteriorRegionIndex() == exteriorRegionID)) {
// Polygon polygon = surface.getPolygons(k);
for (Polygon polygon : surface) {
if (polygonMembaneElementMap != null) {
membraneIndex = polygonMembaneElementMap.get(polygon).getMembraneIndex();
}
Node[] nodes = polygon.getNodes();
if (dimension == 2) {
// ignore z
Vect3d unitNormal = new Vect3d();
polygon.getUnitNormal(unitNormal);
unitNormal.set(unitNormal.getX(), unitNormal.getY(), 0);
int point0 = 0;
Vect3d v0 = new Vect3d(nodes[point0].getX(), nodes[point0].getY(), 0);
int point1 = 1;
Vect3d v1 = null;
for (point1 = 1; point1 < nodes.length; point1++) {
if (v0.getX() != nodes[point1].getX() || v0.getY() != nodes[point1].getY()) {
v1 = new Vect3d(nodes[point1].getX(), nodes[point1].getY(), 0);
break;
}
}
if (v1 == null) {
throw new RuntimeException("failed to generate surface");
}
Vect3d v01 = Vect3d.sub(v1, v0);
Vect3d unit01n = v01.cross(unitNormal);
unit01n.unit();
if (Math.abs(unit01n.getZ() - 1.0) < 1e-6) {
// v0 to v1 opposes vcell surface normal. it's already flipped.
Triangle triangle;
if (surface.getInteriorRegionIndex() == interiorRegionID) {
// we have to flipped it back
triangle = new Triangle(nodes[point1], nodes[point0], null);
} else {
triangle = new Triangle(nodes[point0], nodes[point1], null);
}
triList.add(new TrianglePanel(triLocalCount++, triangleGlobalCount++, membraneIndex, triangle));
} else if (Math.abs(unit01n.getZ() + 1.0) < 1e-6) {
// v0 to v1 is in direction of vcell surface normal.
Triangle triangle;
if (surface.getInteriorRegionIndex() == interiorRegionID) {
triangle = new Triangle(nodes[point0], nodes[point1], null);
} else {
triangle = new Triangle(nodes[point1], nodes[point0], null);
}
triList.add(new TrianglePanel(triLocalCount++, triangleGlobalCount++, membraneIndex, triangle));
} else {
throw new RuntimeException("failed to generate surface");
}
} else if (dimension == 3) {
Triangle triangle1;
Triangle triangle2;
if (surface.getInteriorRegionIndex() == interiorRegionID) {
// interior
triangle1 = new Triangle(nodes[0], nodes[1], nodes[2]);
triangle2 = new Triangle(nodes[0], nodes[2], nodes[3]);
} else {
triangle1 = new Triangle(nodes[2], nodes[1], nodes[0]);
triangle2 = new Triangle(nodes[3], nodes[2], nodes[0]);
}
triList.add(new TrianglePanel(triLocalCount++, triangleGlobalCount++, membraneIndex, triangle1));
triList.add(new TrianglePanel(triLocalCount++, triangleGlobalCount++, membraneIndex, triangle2));
}
}
}
}
}
// add triangles to node hash
for (TrianglePanel triPanel : triList) {
for (Node node : triPanel.triangle.getNodes()) {
if (node == null) {
continue;
}
Set<String> triNameSet = nodeTriMap.get(node);
if (triNameSet == null) {
triNameSet = new HashSet<String>();
nodeTriMap.put(node, triNameSet);
}
triNameSet.add(triPanel.name);
}
}
SubVolume[] adjacentSubvolums = surfaceClass.getAdjacentSubvolumes().toArray(new SubVolume[0]);
CompartmentSubDomain csd0 = simulation.getMathDescription().getCompartmentSubDomain(adjacentSubvolums[0].getName());
CompartmentSubDomain csd1 = simulation.getMathDescription().getCompartmentSubDomain(adjacentSubvolums[1].getName());
MembraneSubDomain membraneSubDomain = simulation.getMathDescription().getMembraneSubDomain(csd0, csd1);
membraneSubdomainTriangleMap.put(membraneSubDomain, triList);
final boolean initialMoleculesOnMembrane = (closestTriangles != null);
if (initialMoleculesOnMembrane) {
findClosestTriangles(membraneSubDomain, triList);
}
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.start_surface + " " + surfaceClass.getName());
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + "(" + SmoldynVCellMapper.SmoldynKeyword.all + ") " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.reflect);
// printWriter.println(SmoldynKeyword.action + " " + SmoldynKeyword.all + "(" + SmoldynKeyword.up + ") " + SmoldynKeyword.both + " " + SmoldynKeyword.reflect);
Color c = colorForSurface(sci);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.color + " " + SmoldynVCellMapper.SmoldynKeyword.both + " " + c.getRed() / 255.0 + " " + c.getGreen() / 255.0 + " " + c.getBlue() / 255.0 + " 0.1");
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.front + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.back + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_panels + " " + SmoldynVCellMapper.SmoldynKeyword.tri + " " + triList.size());
for (TrianglePanel trianglePanel : triList) {
Triangle triangle = trianglePanel.triangle;
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.tri);
switch(dimension) {
case 1:
printWriter.print(" " + triangle.getNodes(0).getX());
break;
case 2:
printWriter.print(" " + triangle.getNodes(0).getX() + " " + triangle.getNodes(0).getY());
printWriter.print(" " + triangle.getNodes(1).getX() + " " + triangle.getNodes(1).getY());
break;
case 3:
for (Node node : triangle.getNodes()) {
printWriter.print(" " + node.getX() + " " + node.getY() + " " + node.getZ());
}
break;
}
printWriter.println(" " + trianglePanel.name);
}
for (TrianglePanel triPanel : triList) {
neighborsForCurrentNode.clear();
for (Node node : triPanel.triangle.getNodes()) {
if (node == null) {
continue;
}
neighborsForCurrentNode.addAll(nodeTriMap.get(node));
}
neighborsForCurrentNode.remove(triPanel.name);
// printWriter.print(SmoldynKeyword.neighbors + " " +triPanel.name);
// to allow smoldyn read line length as 256, chop the neighbors to multiple lines
int maxNeighborCount = 4;
//
int count = 0;
for (String neigh : neighborsForCurrentNode) {
if (count % maxNeighborCount == 0) {
printWriter.println();
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.neighbors + " " + triPanel.name);
}
printWriter.print(" " + neigh);
count++;
}
}
printWriter.println();
printWriter.println(SmoldynVCellMapper.SmoldynKeyword.end_surface);
printWriter.println();
}
}
}
Also used :
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
SortedSet(java.util.SortedSet)
TreeSet(java.util.TreeSet)
HashSet(java.util.HashSet)
Set(java.util.Set)
GeometrySurfaceDescription(cbit.vcell.geometry.surface.GeometrySurfaceDescription)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
HashMap(java.util.HashMap)
Node(cbit.vcell.geometry.surface.Node)
ArrayList(java.util.ArrayList)
Triangle(cbit.vcell.geometry.surface.Triangle)
Surface(cbit.vcell.geometry.surface.Surface)
GeometrySpec(cbit.vcell.geometry.GeometrySpec)
SubVolume(cbit.vcell.geometry.SubVolume)
TreeSet(java.util.TreeSet)
Polygon(cbit.vcell.geometry.surface.Polygon)
SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion)
SurfaceCollection(cbit.vcell.geometry.surface.SurfaceCollection)
Color(java.awt.Color)
VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion)
VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion)
SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion)
GeometricRegion(cbit.vcell.geometry.surface.GeometricRegion)
Vect3d(cbit.vcell.render.Vect3d)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SolverException(cbit.vcell.solver.SolverException)
Example 4 with Node
use of cbit.vcell.geometry.surface.Node in project vcell by virtualcell.
the class RasterExporter method calculateCommonNeighbor.
private static CommonNeighb calculateCommonNeighbor(CartesianMesh mesh, RegionImage regionImage, MembraneElement membraneElementN0, MembraneElement membraneElementN1, int parentMembrIndex) {
// find common node
Quadrilateral quadp = (Quadrilateral) regionImage.getSurfacecollection().getSurfaces(regionImage.getQuadIndexToSurfAndFace().get(parentMembrIndex).getSurf()).getPolygons(regionImage.getQuadIndexToSurfAndFace().get(parentMembrIndex).getFace());
Quadrilateral quad0 = (Quadrilateral) regionImage.getSurfacecollection().getSurfaces(regionImage.getQuadIndexToSurfAndFace().get(membraneElementN0.getMembraneIndex()).getSurf()).getPolygons(regionImage.getQuadIndexToSurfAndFace().get(membraneElementN0.getMembraneIndex()).getFace());
Quadrilateral quad1 = (Quadrilateral) regionImage.getSurfacecollection().getSurfaces(regionImage.getQuadIndexToSurfAndFace().get(membraneElementN1.getMembraneIndex()).getSurf()).getPolygons(regionImage.getQuadIndexToSurfAndFace().get(membraneElementN1.getMembraneIndex()).getFace());
List<Node> nodesp = Arrays.asList(quadp.getNodes());
List<Node> nodes0 = Arrays.asList(quad0.getNodes());
List<Node> nodes1 = Arrays.asList(quad1.getNodes());
List<Node> intersectN = new ArrayList<>(nodes0);
intersectN.retainAll(nodes1);
intersectN.retainAll(nodesp);
if (intersectN.size() != 1) {
System.out.println("Couldn't find parent node");
return null;
}
ArrayList<Integer> commonNeighbors = new ArrayList<>();
findNeighborNeighborsWithNode(regionImage, membraneElementN0, intersectN.get(0), parentMembrIndex, commonNeighbors);
findNeighborNeighborsWithNode(regionImage, membraneElementN1, intersectN.get(0), parentMembrIndex, commonNeighbors);
return new CommonNeighb(intersectN.get(0), commonNeighbors);
}
Also used :
Quadrilateral(cbit.vcell.geometry.surface.Quadrilateral)
Node(cbit.vcell.geometry.surface.Node)
ArrayList(java.util.ArrayList)
Example 5 with Node
use of cbit.vcell.geometry.surface.Node in project vcell by virtualcell.
the class GeometryFileWriter method write.
/**
* Insert the method's description here.
* Creation date: (7/19/2004 10:54:30 AM)
* @param geometrySurfaceDescription cbit.vcell.geometry.surface.GeometrySurfaceDescription
* @throws IOException
*/
public static void write(Writer writer, Geometry resampledGeometry) throws IOException {
//
// "name" name
// "dimension" dimension
// "extent" extentx extenty extentz
// "origin" originx originy originz
// "volumeRegions" num
// name totalVolume featureHandle
// "membraneRegions" num
// name totalArea volumeRegionIndex1 volumeRegionIndex2
// "volumeSamples" numX, numY, numZ
// uncompressed regionIndexs for each volume element
// compressed regionIndexs for each volume element
// "nodes" num
// nodeIndex x y z
// "cells" num
// cellIndex patchIndex node1 node2 node3 node4
// "celldata"
// insideVolumeIndex outsideVolumeIndex area normalx normaly normalz
//
//
// When we are writing volume regions, we sort regions so that ID is equal to index
//
writer.write("name " + resampledGeometry.getName() + "\n");
writer.write("dimension " + resampledGeometry.getDimension() + "\n");
org.vcell.util.Extent extent = resampledGeometry.getExtent();
org.vcell.util.Origin origin = resampledGeometry.getOrigin();
switch(resampledGeometry.getDimension()) {
case 1:
writer.write("size " + extent.getX() + "\n");
writer.write("origin " + origin.getX() + "\n");
break;
case 2:
writer.write("size " + extent.getX() + " " + extent.getY() + "\n");
writer.write("origin " + origin.getX() + " " + origin.getY() + "\n");
break;
case 3:
writer.write("size " + extent.getX() + " " + extent.getY() + " " + extent.getZ() + "\n");
writer.write("origin " + origin.getX() + " " + origin.getY() + " " + origin.getZ() + "\n");
break;
}
GeometrySurfaceDescription geoSurfaceDesc = resampledGeometry.getGeometrySurfaceDescription();
RegionImage regionImage = geoSurfaceDesc.getRegionImage();
SurfaceCollection surfaceCollection = geoSurfaceDesc.getSurfaceCollection();
GeometricRegion[] geometricRegions = geoSurfaceDesc.getGeometricRegions();
int numVolumeRegions = 0;
int numMembraneRegions = 0;
Vector<VolumeGeometricRegion> volRegionList = new Vector<VolumeGeometricRegion>();
if (geometricRegions != null) {
for (int i = 0; i < geometricRegions.length; i++) {
if (geometricRegions[i] instanceof VolumeGeometricRegion) {
numVolumeRegions++;
volRegionList.add((VolumeGeometricRegion) geometricRegions[i]);
} else if (geometricRegions[i] instanceof SurfaceGeometricRegion) {
numMembraneRegions++;
}
}
}
//
// get ordered array of volume regions (where "id" == index into array)... fail if impossible
//
java.util.Collections.sort(volRegionList, new Comparator<VolumeGeometricRegion>() {
public int compare(VolumeGeometricRegion reg1, VolumeGeometricRegion reg2) {
if (reg1.getRegionID() < reg2.getRegionID()) {
return -1;
} else if (reg1.getRegionID() > reg2.getRegionID()) {
return 1;
} else {
return 0;
}
}
public boolean equals(Object obj) {
return this == obj;
}
});
VolumeGeometricRegion[] volRegions = (VolumeGeometricRegion[]) org.vcell.util.BeanUtils.getArray(volRegionList, VolumeGeometricRegion.class);
writer.write("volumeRegions " + numVolumeRegions + "\n");
for (int i = 0; i < volRegions.length; i++) {
if (volRegions[i].getRegionID() != i) {
throw new RuntimeException("Region ID != Region Index, they must be the same!");
}
writer.write(volRegions[i].getName() + " " + volRegions[i].getSize() + " " + volRegions[i].getSubVolume().getHandle() + "\n");
}
writer.write("membraneRegions " + numMembraneRegions + "\n");
if (geometricRegions != null) {
for (int i = 0; i < geometricRegions.length; i++) {
if (geometricRegions[i] instanceof SurfaceGeometricRegion) {
SurfaceGeometricRegion surfaceRegion = (SurfaceGeometricRegion) geometricRegions[i];
GeometricRegion[] neighbors = surfaceRegion.getAdjacentGeometricRegions();
VolumeGeometricRegion insideRegion = (VolumeGeometricRegion) neighbors[0];
VolumeGeometricRegion outsideRegion = (VolumeGeometricRegion) neighbors[1];
writer.write(surfaceRegion.getName() + " " + surfaceRegion.getSize() + " " + insideRegion.getRegionID() + " " + outsideRegion.getRegionID() + "\n");
}
}
}
//
// write volume samples
//
ISize volumeSampleSize = geoSurfaceDesc.getVolumeSampleSize();
switch(resampledGeometry.getDimension()) {
case 1:
writer.write("volumeSamples " + volumeSampleSize.getX() + "\n");
break;
case 2:
writer.write("volumeSamples " + volumeSampleSize.getX() + " " + volumeSampleSize.getY() + "\n");
break;
case 3:
writer.write("volumeSamples " + volumeSampleSize.getX() + " " + volumeSampleSize.getY() + " " + volumeSampleSize.getZ() + "\n");
break;
}
// regionImage
if (regionImage != null) {
if (regionImage.getNumRegions() > 65536) {
throw new RuntimeException("cannot process a geometry with more than 65536 volume regions");
}
byte[] uncompressedRegionIDs = new byte[2 * regionImage.getNumX() * regionImage.getNumY() * regionImage.getNumZ()];
for (int i = 0, j = 0; i < uncompressedRegionIDs.length; i += 2, j++) {
int regindex = regionImage.getRegionInfoFromOffset(j).getRegionIndex();
uncompressedRegionIDs[i] = (byte) (regindex & 0x000000ff);
uncompressedRegionIDs[i + 1] = (byte) ((regindex & 0x0000ff00) >> 8);
}
ByteArrayOutputStream bos = new ByteArrayOutputStream();
DeflaterOutputStream dos = new DeflaterOutputStream(bos);
dos.write(uncompressedRegionIDs, 0, uncompressedRegionIDs.length);
dos.close();
byte[] compressedRegionIDs = bos.toByteArray();
writer.write(org.vcell.util.Hex.toString(compressedRegionIDs) + "\n");
} else {
writer.write("\n");
}
//
if (surfaceCollection == null) {
throw new RuntimeException("geometry is not updated");
}
int numCells = surfaceCollection.getTotalPolygonCount();
writer.write("cells " + numCells + "\n");
// "celldata"
// insideVolumeIndex outsideVolumeIndex area normalx normaly normalz
//
int cellID = 0;
int dimension = resampledGeometry.getDimension();
double correctCoeff = 1;
if (dimension == 1) {
correctCoeff = extent.getY() * extent.getZ();
} else if (dimension == 2) {
correctCoeff = extent.getZ();
}
if (surfaceCollection != null) {
for (int i = 0; i < surfaceCollection.getSurfaceCount(); i++) {
Surface surface = surfaceCollection.getSurfaces(i);
int region1Outside = 0;
int region1Inside = 0;
for (int j = 0; j < surface.getPolygonCount(); j++) {
Quadrilateral polygon = (Quadrilateral) surface.getPolygons(j);
Node[] node = polygon.getNodes();
cbit.vcell.render.Vect3d elementCoord = new cbit.vcell.render.Vect3d();
int nodesOnBoundary = 0;
for (int k = 0; k < node.length; k++) {
if (!node[k].getMoveX() || (dimension > 1 && !node[k].getMoveY()) || (dimension == 3 && !node[k].getMoveZ())) {
nodesOnBoundary++;
}
}
if (nodesOnBoundary == 0) {
for (int k = 0; k < node.length; k++) {
elementCoord.add(new cbit.vcell.render.Vect3d(node[k].getX(), node[k].getY(), node[k].getZ()));
}
elementCoord.scale(0.25);
} else if (nodesOnBoundary == 2) {
for (int k = 0; k < node.length; k++) {
if (!node[k].getMoveX() || !node[k].getMoveY() || !node[k].getMoveZ()) {
elementCoord.add(new cbit.vcell.render.Vect3d(node[k].getX(), node[k].getY(), node[k].getZ()));
}
}
elementCoord.scale(0.5);
} else if (nodesOnBoundary == 3) {
for (int k = 0; k < node.length; k++) {
if (!node[k].getMoveX() && !node[k].getMoveY() || !node[k].getMoveY() && !node[k].getMoveZ() || !node[k].getMoveX() && !node[k].getMoveZ()) {
elementCoord.set(node[k].getX(), node[k].getY(), node[k].getZ());
}
}
} else {
throw new RuntimeException("Unexcepted number of nodes on boundary for a polygon: " + nodesOnBoundary);
}
cbit.vcell.render.Vect3d unitNormal = new cbit.vcell.render.Vect3d();
polygon.getUnitNormal(unitNormal);
int volNeighbor1Region = regionImage.getRegionInfoFromOffset(polygon.getVolIndexNeighbor1()).getRegionIndex();
int volNeighbor2Region = regionImage.getRegionInfoFromOffset(polygon.getVolIndexNeighbor2()).getRegionIndex();
if (surface.getExteriorRegionIndex() == volNeighbor1Region && surface.getInteriorRegionIndex() == volNeighbor2Region) {
region1Outside++;
}
if (surface.getExteriorRegionIndex() == volNeighbor2Region && surface.getInteriorRegionIndex() == volNeighbor1Region) {
region1Inside++;
}
writer.write(cellID + " " + polygon.getVolIndexNeighbor1() + " " + polygon.getVolIndexNeighbor2() + " " + polygon.getArea() / correctCoeff + " " + elementCoord.getX() + " " + elementCoord.getY() + " " + elementCoord.getZ() + " " + unitNormal.getX() + " " + unitNormal.getY() + " " + unitNormal.getZ() + "\n");
cellID++;
}
if (region1Inside != surface.getPolygonCount() && region1Outside != surface.getPolygonCount()) {
throw new RuntimeException("Volume neighbor regions not consistent: [total, inside, outside]=" + surface.getPolygonCount() + "," + region1Inside + "," + region1Outside + "]");
}
}
}
}
Also used :
GeometrySurfaceDescription(cbit.vcell.geometry.surface.GeometrySurfaceDescription)
ISize(org.vcell.util.ISize)
Node(cbit.vcell.geometry.surface.Node)
Surface(cbit.vcell.geometry.surface.Surface)
Quadrilateral(cbit.vcell.geometry.surface.Quadrilateral)
DeflaterOutputStream(java.util.zip.DeflaterOutputStream)
Vector(java.util.Vector)
SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion)
SurfaceCollection(cbit.vcell.geometry.surface.SurfaceCollection)
ByteArrayOutputStream(java.io.ByteArrayOutputStream)
VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion)
VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion)
SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion)
GeometricRegion(cbit.vcell.geometry.surface.GeometricRegion)
RegionImage(cbit.vcell.geometry.RegionImage)
Aggregations
Node (cbit.vcell.geometry.surface.Node)9
Surface (cbit.vcell.geometry.surface.Surface)5
SurfaceCollection (cbit.vcell.geometry.surface.SurfaceCollection)5
Polygon (cbit.vcell.geometry.surface.Polygon)4
Quadrilateral (cbit.vcell.geometry.surface.Quadrilateral)4
ArrayList (java.util.ArrayList)4
GeometricRegion (cbit.vcell.geometry.surface.GeometricRegion)3
GeometrySurfaceDescription (cbit.vcell.geometry.surface.GeometrySurfaceDescription)3
SurfaceGeometricRegion (cbit.vcell.geometry.surface.SurfaceGeometricRegion)3
VolumeGeometricRegion (cbit.vcell.geometry.surface.VolumeGeometricRegion)3
TreeSet (java.util.TreeSet)3
GeometrySpec (cbit.vcell.geometry.GeometrySpec)2
RegionImage (cbit.vcell.geometry.RegionImage)2
SubVolume (cbit.vcell.geometry.SubVolume)2
SurfaceClass (cbit.vcell.geometry.SurfaceClass)2
Triangle (cbit.vcell.geometry.surface.Triangle)2
CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)2
MembraneSubDomain (cbit.vcell.math.MembraneSubDomain)2
Vect3d (cbit.vcell.render.Vect3d)2
SolverException (cbit.vcell.solver.SolverException)2
