Examples with GeometrySpec cbit.vcell.geometry.GeometrySpec used on opensource projects

Example 1 with GeometrySpec

use of cbit.vcell.geometry.GeometrySpec 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();
    }
}

Example 2 with GeometrySpec

use of cbit.vcell.geometry.GeometrySpec 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();
        }
    }
}

Example 3 with GeometrySpec

use of cbit.vcell.geometry.GeometrySpec in project vcell by virtualcell.

the class SBMLImporter method addGeometry.

protected void addGeometry() {
    // get a Geometry object via SpatialModelPlugin object.
    org.sbml.jsbml.ext.spatial.Geometry sbmlGeometry = getSbmlGeometry();
    if (sbmlGeometry == null) {
        return;
    }
    int dimension = 0;
    Origin vcOrigin = null;
    Extent vcExtent = null;
    {
        // local code block
        // get a CoordComponent object via the Geometry object.
        ListOf<CoordinateComponent> listOfCoordComps = sbmlGeometry.getListOfCoordinateComponents();
        if (listOfCoordComps == null) {
            throw new RuntimeException("Cannot have 0 coordinate compartments in geometry");
        }
        // coord component
        double ox = 0.0;
        double oy = 0.0;
        double oz = 0.0;
        double ex = 1.0;
        double ey = 1.0;
        double ez = 1.0;
        for (CoordinateComponent coordComponent : listOfCoordComps) {
            double minValue = coordComponent.getBoundaryMinimum().getValue();
            double maxValue = coordComponent.getBoundaryMaximum().getValue();
            switch(coordComponent.getType()) {
                case cartesianX:
                    {
                        ox = minValue;
                        ex = maxValue - minValue;
                        break;
                    }
                case cartesianY:
                    {
                        oy = minValue;
                        ey = maxValue - minValue;
                        break;
                    }
                case cartesianZ:
                    {
                        oz = minValue;
                        ez = maxValue - minValue;
                        break;
                    }
            }
            dimension++;
        }
        vcOrigin = new Origin(ox, oy, oz);
        vcExtent = new Extent(ex, ey, ez);
    }
    // from geometry definition, find out which type of geometry : image or
    // analytic or CSG
    AnalyticGeometry analyticGeometryDefinition = null;
    CSGeometry csGeometry = null;
    SampledFieldGeometry segmentedSampledFieldGeometry = null;
    SampledFieldGeometry distanceMapSampledFieldGeometry = null;
    ParametricGeometry parametricGeometry = null;
    for (int i = 0; i < sbmlGeometry.getListOfGeometryDefinitions().size(); i++) {
        GeometryDefinition gd_temp = sbmlGeometry.getListOfGeometryDefinitions().get(i);
        if (!gd_temp.isSetIsActive()) {
            continue;
        }
        if (gd_temp instanceof AnalyticGeometry) {
            analyticGeometryDefinition = (AnalyticGeometry) gd_temp;
        } else if (gd_temp instanceof SampledFieldGeometry) {
            SampledFieldGeometry sfg = (SampledFieldGeometry) gd_temp;
            String sfn = sfg.getSampledField();
            ListOf<SampledField> sampledFields = sbmlGeometry.getListOfSampledFields();
            if (sampledFields.size() > 1) {
                throw new RuntimeException("only one sampled field supported");
            }
            InterpolationKind ik = sampledFields.get(0).getInterpolationType();
            switch(ik) {
                case linear:
                    distanceMapSampledFieldGeometry = sfg;
                    break;
                case nearestneighbor:
                    segmentedSampledFieldGeometry = sfg;
                    break;
                default:
                    lg.warn("Unsupported " + sampledFields.get(0).getName() + " interpolation type " + ik);
            }
        } else if (gd_temp instanceof CSGeometry) {
            csGeometry = (CSGeometry) gd_temp;
        } else if (gd_temp instanceof ParametricGeometry) {
            parametricGeometry = (ParametricGeometry) gd_temp;
        } else {
            throw new RuntimeException("unsupported geometry definition type " + gd_temp.getClass().getSimpleName());
        }
    }
    if (analyticGeometryDefinition == null && segmentedSampledFieldGeometry == null && distanceMapSampledFieldGeometry == null && csGeometry == null) {
        throw new SBMLImportException("VCell supports only Analytic, Image based (segmentd or distance map) or Constructed Solid Geometry at this time.");
    }
    GeometryDefinition selectedGeometryDefinition = null;
    if (csGeometry != null) {
        selectedGeometryDefinition = csGeometry;
    } else if (analyticGeometryDefinition != null) {
        selectedGeometryDefinition = analyticGeometryDefinition;
    } else if (segmentedSampledFieldGeometry != null) {
        selectedGeometryDefinition = segmentedSampledFieldGeometry;
    } else if (distanceMapSampledFieldGeometry != null) {
        selectedGeometryDefinition = distanceMapSampledFieldGeometry;
    } else if (parametricGeometry != null) {
        selectedGeometryDefinition = parametricGeometry;
    } else {
        throw new SBMLImportException("no geometry definition found");
    }
    Geometry vcGeometry = null;
    if (selectedGeometryDefinition == analyticGeometryDefinition || selectedGeometryDefinition == csGeometry) {
        vcGeometry = new Geometry("spatialGeom", dimension);
    } else if (selectedGeometryDefinition == distanceMapSampledFieldGeometry || selectedGeometryDefinition == segmentedSampledFieldGeometry) {
        SampledFieldGeometry sfg = (SampledFieldGeometry) selectedGeometryDefinition;
        // get image from sampledFieldGeometry
        // get a sampledVol object via the listOfSampledVol (from
        // SampledGeometry) object.
        // gcw gcw gcw
        String sfn = sfg.getSampledField();
        SampledField sf = null;
        for (SampledField sampledField : sbmlGeometry.getListOfSampledFields()) {
            if (sampledField.getSpatialId().equals(sfn)) {
                sf = sampledField;
            }
        }
        int numX = sf.getNumSamples1();
        int numY = sf.getNumSamples2();
        int numZ = sf.getNumSamples3();
        int[] samples = new int[sf.getSamplesLength()];
        StringTokenizer tokens = new StringTokenizer(sf.getSamples(), " ");
        int count = 0;
        while (tokens.hasMoreTokens()) {
            int sample = Integer.parseInt(tokens.nextToken());
            samples[count++] = sample;
        }
        byte[] imageInBytes = new byte[samples.length];
        if (selectedGeometryDefinition == distanceMapSampledFieldGeometry) {
            // 
            for (int i = 0; i < imageInBytes.length; i++) {
                // if (interpolation(samples[i])<0){
                if (samples[i] < 0) {
                    imageInBytes[i] = -1;
                } else {
                    imageInBytes[i] = 1;
                }
            }
        } else {
            for (int i = 0; i < imageInBytes.length; i++) {
                imageInBytes[i] = (byte) samples[i];
            }
        }
        try {
            // System.out.println("ident " + sf.getId() + " " + sf.getName());
            VCImage vcImage = null;
            CompressionKind ck = sf.getCompression();
            DataKind dk = sf.getDataType();
            if (ck == CompressionKind.deflated) {
                vcImage = new VCImageCompressed(null, imageInBytes, vcExtent, numX, numY, numZ);
            } else {
                switch(dk) {
                    case UINT8:
                    case UINT16:
                    case UINT32:
                        vcImage = new VCImageUncompressed(null, imageInBytes, vcExtent, numX, numY, numZ);
                    default:
                }
            }
            if (vcImage == null) {
                throw new SbmlException("Unsupported type combination " + ck + ", " + dk + " for sampled field " + sf.getName());
            }
            vcImage.setName(sf.getId());
            ListOf<SampledVolume> sampledVolumes = sfg.getListOfSampledVolumes();
            final int numSampledVols = sampledVolumes.size();
            if (numSampledVols == 0) {
                throw new RuntimeException("Cannot have 0 sampled volumes in sampledField (image_based) geometry");
            }
            // check to see if values are uniquely integer , add set up scaling if necessary
            double scaleFactor = checkPixelScaling(sampledVolumes, 1);
            if (scaleFactor != 1) {
                double checkScaleFactor = checkPixelScaling(sampledVolumes, scaleFactor);
                VCAssert.assertTrue(checkScaleFactor != scaleFactor, "Scale factor check failed");
            }
            VCPixelClass[] vcpixelClasses = new VCPixelClass[numSampledVols];
            // get pixel classes for geometry
            for (int i = 0; i < numSampledVols; i++) {
                SampledVolume sVol = sampledVolumes.get(i);
                // from subVolume, get pixelClass?
                final int scaled = (int) (scaleFactor * sVol.getSampledValue());
                vcpixelClasses[i] = new VCPixelClass(null, sVol.getDomainType(), scaled);
            }
            vcImage.setPixelClasses(vcpixelClasses);
            // now create image geometry
            vcGeometry = new Geometry("spatialGeom", vcImage);
        } catch (Exception e) {
            e.printStackTrace(System.out);
            throw new RuntimeException("Unable to create image from SampledFieldGeometry : " + e.getMessage());
        }
    }
    GeometrySpec vcGeometrySpec = vcGeometry.getGeometrySpec();
    vcGeometrySpec.setOrigin(vcOrigin);
    try {
        vcGeometrySpec.setExtent(vcExtent);
    } catch (PropertyVetoException e) {
        e.printStackTrace(System.out);
        throw new SBMLImportException("Unable to set extent on VC geometry : " + e.getMessage(), e);
    }
    // get listOfDomainTypes via the Geometry object.
    ListOf<DomainType> listOfDomainTypes = sbmlGeometry.getListOfDomainTypes();
    if (listOfDomainTypes == null || listOfDomainTypes.size() < 1) {
        throw new SBMLImportException("Cannot have 0 domainTypes in geometry");
    }
    // get a listOfDomains via the Geometry object.
    ListOf<Domain> listOfDomains = sbmlGeometry.getListOfDomains();
    if (listOfDomains == null || listOfDomains.size() < 1) {
        throw new SBMLImportException("Cannot have 0 domains in geometry");
    }
    // ListOfGeometryDefinitions listOfGeomDefns =
    // sbmlGeometry.getListOfGeometryDefinitions();
    // if ((listOfGeomDefns == null) ||
    // (sbmlGeometry.getNumGeometryDefinitions() > 1)) {
    // throw new
    // RuntimeException("Can have only 1 geometry definition in geometry");
    // }
    // use the boolean bAnalytic to create the right kind of subvolume.
    // First match the somVol=domainTypes for spDim=3. Deal witl spDim=2
    // afterwards.
    GeometrySurfaceDescription vcGsd = vcGeometry.getGeometrySurfaceDescription();
    Vector<DomainType> surfaceClassDomainTypesVector = new Vector<DomainType>();
    try {
        for (DomainType dt : listOfDomainTypes) {
            if (dt.getSpatialDimensions() == 3) {
                // subvolume
                if (selectedGeometryDefinition == analyticGeometryDefinition) {
                    // will set expression later - when reading in Analytic
                    // Volumes in GeometryDefinition
                    vcGeometrySpec.addSubVolume(new AnalyticSubVolume(dt.getId(), new Expression(1.0)));
                } else {
                // add SubVolumes later for CSG and Image-based
                }
            } else if (dt.getSpatialDimensions() == 2) {
                surfaceClassDomainTypesVector.add(dt);
            }
        }
        // analytic vol is needed to get the expression for subVols
        if (selectedGeometryDefinition == analyticGeometryDefinition) {
            // get an analyticVol object via the listOfAnalyticVol (from
            // AnalyticGeometry) object.
            ListOf<AnalyticVolume> aVolumes = analyticGeometryDefinition.getListOfAnalyticVolumes();
            if (aVolumes.size() < 1) {
                throw new SBMLImportException("Cannot have 0 Analytic volumes in analytic geometry");
            }
            for (AnalyticVolume analyticVol : aVolumes) {
                // get subVol from VC geometry using analyticVol spatialId;
                // set its expr using analyticVol's math.
                SubVolume vcSubvolume = vcGeometrySpec.getSubVolume(analyticVol.getDomainType());
                CastInfo<AnalyticSubVolume> ci = BeanUtils.attemptCast(AnalyticSubVolume.class, vcSubvolume);
                if (!ci.isGood()) {
                    throw new RuntimeException("analytic volume '" + analyticVol.getId() + "' does not map to any VC subvolume.");
                }
                AnalyticSubVolume asv = ci.get();
                try {
                    Expression subVolExpr = getExpressionFromFormula(analyticVol.getMath());
                    asv.setExpression(subVolExpr);
                } catch (ExpressionException e) {
                    e.printStackTrace(System.out);
                    throw new SBMLImportException("Unable to set expression on subVolume '" + asv.getName() + "'. " + e.getMessage(), e);
                }
            }
        }
        SampledFieldGeometry sfg = BeanUtils.downcast(SampledFieldGeometry.class, selectedGeometryDefinition);
        if (sfg != null) {
            ListOf<SampledVolume> sampledVolumes = sfg.getListOfSampledVolumes();
            int numSampledVols = sampledVolumes.size();
            if (numSampledVols == 0) {
                throw new SBMLImportException("Cannot have 0 sampled volumes in sampledField (image_based) geometry");
            }
            VCPixelClass[] vcpixelClasses = new VCPixelClass[numSampledVols];
            ImageSubVolume[] vcImageSubVols = new ImageSubVolume[numSampledVols];
            // get pixel classes for geometry
            int idx = 0;
            for (SampledVolume sVol : sampledVolumes) {
                // from subVolume, get pixelClass?
                final String name = sVol.getDomainType();
                final int pixelValue = SBMLUtils.ignoreZeroFraction(sVol.getSampledValue());
                VCPixelClass pc = new VCPixelClass(null, name, pixelValue);
                vcpixelClasses[idx] = pc;
                // Create the new Image SubVolume - use index of this for
                // loop as 'handle' for ImageSubVol?
                ImageSubVolume isv = new ImageSubVolume(null, pc, idx);
                isv.setName(name);
                vcImageSubVols[idx++] = isv;
            }
            vcGeometry.getGeometrySpec().setSubVolumes(vcImageSubVols);
        }
        if (selectedGeometryDefinition == csGeometry) {
            ListOf<org.sbml.jsbml.ext.spatial.CSGObject> listOfcsgObjs = csGeometry.getListOfCSGObjects();
            ArrayList<org.sbml.jsbml.ext.spatial.CSGObject> sbmlCSGs = new ArrayList<org.sbml.jsbml.ext.spatial.CSGObject>(listOfcsgObjs);
            // we want the CSGObj with highest ordinal to be the first
            // element in the CSG subvols array.
            Collections.sort(sbmlCSGs, new Comparator<org.sbml.jsbml.ext.spatial.CSGObject>() {

                @Override
                public int compare(org.sbml.jsbml.ext.spatial.CSGObject lhs, org.sbml.jsbml.ext.spatial.CSGObject rhs) {
                    // minus one to reverse sort
                    return -1 * Integer.compare(lhs.getOrdinal(), rhs.getOrdinal());
                }
            });
            int n = sbmlCSGs.size();
            CSGObject[] vcCSGSubVolumes = new CSGObject[n];
            for (int i = 0; i < n; i++) {
                org.sbml.jsbml.ext.spatial.CSGObject sbmlCSGObject = sbmlCSGs.get(i);
                CSGObject vcellCSGObject = new CSGObject(null, sbmlCSGObject.getDomainType(), i);
                vcellCSGObject.setRoot(getVCellCSGNode(sbmlCSGObject.getCSGNode()));
            }
            vcGeometry.getGeometrySpec().setSubVolumes(vcCSGSubVolumes);
        }
        // Call geom.geomSurfDesc.updateAll() to automatically generate
        // surface classes.
        // vcGsd.updateAll();
        vcGeometry.precomputeAll(new GeometryThumbnailImageFactoryAWT(), true, true);
    } catch (Exception e) {
        e.printStackTrace(System.out);
        throw new SBMLImportException("Unable to create VC subVolumes from SBML domainTypes : " + e.getMessage(), e);
    }
    // should now map each SBML domain to right VC geometric region.
    GeometricRegion[] vcGeomRegions = vcGsd.getGeometricRegions();
    ISize sampleSize = vcGsd.getVolumeSampleSize();
    RegionInfo[] regionInfos = vcGsd.getRegionImage().getRegionInfos();
    int numX = sampleSize.getX();
    int numY = sampleSize.getY();
    int numZ = sampleSize.getZ();
    double ox = vcOrigin.getX();
    double oy = vcOrigin.getY();
    double oz = vcOrigin.getZ();
    for (Domain domain : listOfDomains) {
        String domainType = domain.getDomainType();
        InteriorPoint interiorPt = domain.getListOfInteriorPoints().get(0);
        if (interiorPt == null) {
            DomainType currDomainType = null;
            for (DomainType dt : sbmlGeometry.getListOfDomainTypes()) {
                if (dt.getSpatialId().equals(domainType)) {
                    currDomainType = dt;
                }
            }
            if (currDomainType.getSpatialDimensions() == 2) {
                continue;
            }
        }
        Coordinate sbmlInteriorPtCoord = new Coordinate(interiorPt.getCoord1(), interiorPt.getCoord2(), interiorPt.getCoord3());
        for (int j = 0; j < vcGeomRegions.length; j++) {
            if (vcGeomRegions[j] instanceof VolumeGeometricRegion) {
                int regionID = ((VolumeGeometricRegion) vcGeomRegions[j]).getRegionID();
                for (int k = 0; k < regionInfos.length; k++) {
                    // (using gemoRegion regionID).
                    if (regionInfos[k].getRegionIndex() == regionID) {
                        int volIndx = 0;
                        Coordinate nearestPtCoord = null;
                        double minDistance = Double.MAX_VALUE;
                        // represented by SBML 'domain[i]'.
                        for (int z = 0; z < numZ; z++) {
                            for (int y = 0; y < numY; y++) {
                                for (int x = 0; x < numX; x++) {
                                    if (regionInfos[k].isIndexInRegion(volIndx)) {
                                        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;
                                        // for now, find the shortest dist
                                        // coord. Can refine algo later.
                                        Coordinate vcCoord = new Coordinate(coordX, coordY, coordZ);
                                        double distance = sbmlInteriorPtCoord.distanceTo(vcCoord);
                                        if (distance < minDistance) {
                                            minDistance = distance;
                                            nearestPtCoord = vcCoord;
                                        }
                                    }
                                    volIndx++;
                                }
                            // end - for x
                            }
                        // end - for y
                        }
                        // with domain name
                        if (nearestPtCoord != null) {
                            GeometryClass geomClassSBML = vcGeometry.getGeometryClass(domainType);
                            // we know vcGeometryReg[j] is a VolGeomRegion
                            GeometryClass geomClassVC = ((VolumeGeometricRegion) vcGeomRegions[j]).getSubVolume();
                            if (geomClassSBML.compareEqual(geomClassVC)) {
                                vcGeomRegions[j].setName(domain.getId());
                            }
                        }
                    }
                // end if (regInfoIndx = regId)
                }
            // end - for regInfo
            }
        }
    // end for - vcGeomRegions
    }
    // deal with surfaceClass:spDim2-domainTypes
    for (int i = 0; i < surfaceClassDomainTypesVector.size(); i++) {
        DomainType surfaceClassDomainType = surfaceClassDomainTypesVector.elementAt(i);
        // 'surfaceClassDomainType'
        for (Domain d : listOfDomains) {
            if (d.getDomainType().equals(surfaceClassDomainType.getId())) {
                // get the adjacent domains of this 'surface' domain
                // (surface domain + its 2 adj vol domains)
                Set<Domain> adjacentDomainsSet = getAssociatedAdjacentDomains(sbmlGeometry, d);
                // get the domain types of the adjacent domains in SBML and
                // store the corresponding subVol counterparts from VC for
                // adj vol domains
                Vector<SubVolume> adjacentSubVolumesVector = new Vector<SubVolume>();
                Vector<VolumeGeometricRegion> adjVolGeomRegionsVector = new Vector<VolumeGeometricRegion>();
                Iterator<Domain> iterator = adjacentDomainsSet.iterator();
                while (iterator.hasNext()) {
                    Domain dom = iterator.next();
                    DomainType dt = getBySpatialID(sbmlGeometry.getListOfDomainTypes(), dom.getDomainType());
                    if (dt.getSpatialDimensions() == 3) {
                        // for domain type with sp. dim = 3, get
                        // correspoinding subVol from VC geometry.
                        GeometryClass gc = vcGeometry.getGeometryClass(dt.getId());
                        adjacentSubVolumesVector.add((SubVolume) gc);
                        // store volGeomRegions corresponding to this (vol)
                        // geomClass in adjVolGeomRegionsVector : this
                        // should return ONLY 1 region for subVol.
                        GeometricRegion[] geomRegion = vcGsd.getGeometricRegions(gc);
                        adjVolGeomRegionsVector.add((VolumeGeometricRegion) geomRegion[0]);
                    }
                }
                // there should be only 2 subVols in this vector
                if (adjacentSubVolumesVector.size() != 2) {
                    throw new RuntimeException("Cannot have more or less than 2 subvolumes that are adjacent to surface (membrane) '" + d.getId() + "'");
                }
                // get the surface class with these 2 adj subVols. Set its
                // name to that of 'surfaceClassDomainType'
                SurfaceClass surfacClass = vcGsd.getSurfaceClass(adjacentSubVolumesVector.get(0), adjacentSubVolumesVector.get(1));
                surfacClass.setName(surfaceClassDomainType.getSpatialId());
                // get surfaceGeometricRegion that has adjVolGeomRegions as
                // its adjacent vol geom regions and set its name from
                // domain 'd'
                SurfaceGeometricRegion surfaceGeomRegion = getAssociatedSurfaceGeometricRegion(vcGsd, adjVolGeomRegionsVector);
                if (surfaceGeomRegion != null) {
                    surfaceGeomRegion.setName(d.getId());
                }
            }
        // end if - domain.domainType == surfaceClassDomainType
        }
    // end for - numDomains
    }
    // structureMappings in VC from compartmentMappings in SBML
    try {
        // set geometry first and then set structureMappings?
        vcBioModel.getSimulationContext(0).setGeometry(vcGeometry);
        // update simContextName ...
        vcBioModel.getSimulationContext(0).setName(vcBioModel.getSimulationContext(0).getName() + "_" + vcGeometry.getName());
        Model vcModel = vcBioModel.getSimulationContext(0).getModel();
        ModelUnitSystem vcModelUnitSystem = vcModel.getUnitSystem();
        Vector<StructureMapping> structMappingsVector = new Vector<StructureMapping>();
        SpatialCompartmentPlugin cplugin = null;
        for (int i = 0; i < sbmlModel.getNumCompartments(); i++) {
            Compartment c = sbmlModel.getCompartment(i);
            String cname = c.getName();
            cplugin = (SpatialCompartmentPlugin) c.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
            CompartmentMapping compMapping = cplugin.getCompartmentMapping();
            if (compMapping != null) {
                // final String id = compMapping.getId();
                // final String name = compMapping.getName();
                CastInfo<Structure> ci = SBMLHelper.getTypedStructure(Structure.class, vcModel, cname);
                if (ci.isGood()) {
                    Structure struct = ci.get();
                    String domainType = compMapping.getDomainType();
                    GeometryClass geometryClass = vcGeometry.getGeometryClass(domainType);
                    double unitSize = compMapping.getUnitSize();
                    Feature feat = BeanUtils.downcast(Feature.class, struct);
                    if (feat != null) {
                        FeatureMapping featureMapping = new FeatureMapping(feat, vcBioModel.getSimulationContext(0), vcModelUnitSystem);
                        featureMapping.setGeometryClass(geometryClass);
                        if (geometryClass instanceof SubVolume) {
                            featureMapping.getVolumePerUnitVolumeParameter().setExpression(new Expression(unitSize));
                        } else if (geometryClass instanceof SurfaceClass) {
                            featureMapping.getVolumePerUnitAreaParameter().setExpression(new Expression(unitSize));
                        }
                        structMappingsVector.add(featureMapping);
                    } else if (struct instanceof Membrane) {
                        MembraneMapping membraneMapping = new MembraneMapping((Membrane) struct, vcBioModel.getSimulationContext(0), vcModelUnitSystem);
                        membraneMapping.setGeometryClass(geometryClass);
                        if (geometryClass instanceof SubVolume) {
                            membraneMapping.getAreaPerUnitVolumeParameter().setExpression(new Expression(unitSize));
                        } else if (geometryClass instanceof SurfaceClass) {
                            membraneMapping.getAreaPerUnitAreaParameter().setExpression(new Expression(unitSize));
                        }
                        structMappingsVector.add(membraneMapping);
                    }
                }
            }
        }
        StructureMapping[] structMappings = structMappingsVector.toArray(new StructureMapping[0]);
        vcBioModel.getSimulationContext(0).getGeometryContext().setStructureMappings(structMappings);
        // if type from SBML parameter Boundary Condn is not the same as the
        // boundary type of the
        // structureMapping of structure of paramSpContext, set the boundary
        // condn type of the structureMapping
        // to the value of 'type' from SBML parameter Boundary Condn.
        ListOf<Parameter> listOfGlobalParams = sbmlModel.getListOfParameters();
        for (Parameter sbmlGlobalParam : sbmlModel.getListOfParameters()) {
            SpatialParameterPlugin spplugin = (SpatialParameterPlugin) sbmlGlobalParam.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
            ParameterType paramType = spplugin.getParamType();
            if (!(paramType instanceof BoundaryCondition)) {
                continue;
            }
            BoundaryCondition bCondn = (BoundaryCondition) paramType;
            if (bCondn.isSetVariable()) {
                // get the var of boundaryCondn; find appropriate spContext
                // in vcell;
                SpeciesContext paramSpContext = vcBioModel.getSimulationContext(0).getModel().getSpeciesContext(bCondn.getVariable());
                if (paramSpContext != null) {
                    Structure s = paramSpContext.getStructure();
                    StructureMapping sm = vcBioModel.getSimulationContext(0).getGeometryContext().getStructureMapping(s);
                    if (sm != null) {
                        BoundaryConditionType bct = null;
                        switch(bCondn.getType()) {
                            case Dirichlet:
                                {
                                    bct = BoundaryConditionType.DIRICHLET;
                                    break;
                                }
                            case Neumann:
                                {
                                    bct = BoundaryConditionType.NEUMANN;
                                    break;
                                }
                            case Robin_inwardNormalGradientCoefficient:
                            case Robin_sum:
                            case Robin_valueCoefficient:
                            default:
                                throw new RuntimeException("boundary condition type " + bCondn.getType().name() + " not supported");
                        }
                        for (CoordinateComponent coordComp : getSbmlGeometry().getListOfCoordinateComponents()) {
                            if (bCondn.getSpatialRef().equals(coordComp.getBoundaryMinimum().getSpatialId())) {
                                switch(coordComp.getType()) {
                                    case cartesianX:
                                        {
                                            sm.setBoundaryConditionTypeXm(bct);
                                        }
                                    case cartesianY:
                                        {
                                            sm.setBoundaryConditionTypeYm(bct);
                                        }
                                    case cartesianZ:
                                        {
                                            sm.setBoundaryConditionTypeZm(bct);
                                        }
                                }
                            }
                            if (bCondn.getSpatialRef().equals(coordComp.getBoundaryMaximum().getSpatialId())) {
                                switch(coordComp.getType()) {
                                    case cartesianX:
                                        {
                                            sm.setBoundaryConditionTypeXm(bct);
                                        }
                                    case cartesianY:
                                        {
                                            sm.setBoundaryConditionTypeYm(bct);
                                        }
                                    case cartesianZ:
                                        {
                                            sm.setBoundaryConditionTypeZm(bct);
                                        }
                                }
                            }
                        }
                    } else // sm != null
                    {
                        logger.sendMessage(VCLogger.Priority.MediumPriority, VCLogger.ErrorType.OverallWarning, "No structure " + s.getName() + " requested by species context " + paramSpContext.getName());
                    }
                }
            // end if (paramSpContext != null)
            }
        // end if (bCondn.isSetVar())
        }
        // end for (sbmlModel.numParams)
        vcBioModel.getSimulationContext(0).getGeometryContext().refreshStructureMappings();
        vcBioModel.getSimulationContext(0).refreshSpatialObjects();
    } catch (Exception e) {
        e.printStackTrace(System.out);
        throw new SBMLImportException("Unable to create VC structureMappings from SBML compartment mappings : " + e.getMessage(), e);
    }
}

Example 4 with GeometrySpec

use of cbit.vcell.geometry.GeometrySpec in project vcell by virtualcell.

the class GeometryWindowManager method propertyChange.

/**
 * This method gets called when a bound property is changed.
 * @param evt A PropertyChangeEvent object describing the event source
 *   	and the property that has changed.
 */
public void propertyChange(java.beans.PropertyChangeEvent evt) {
    if (evt.getSource() instanceof GeometrySpec) {
        getGeometryEditor().setToggleButtonSelected("Surface Viewer", false);
        // do this because button.setSelected does not fire action event
        ChildWindow childWindow = ChildWindowManager.findChildWindowManager(getComponent()).getChildWindowFromContentPane(surfaceViewer);
        if (childWindow != null && childWindow.isShowing()) {
            surfaceViewerButtonPressed(false);
        }
    }
}

Example 5 with GeometrySpec

use of cbit.vcell.geometry.GeometrySpec in project vcell by virtualcell.

the class FiniteVolumeFileWriter method writeChomboSpec.

private void writeChomboSpec() throws ExpressionException, SolverException, PropertyVetoException, ClassNotFoundException, IOException, GeometryException, ImageException {
    if (!bChomboSolver) {
        return;
    }
    GeometrySpec geometrySpec = resampledGeometry.getGeometrySpec();
    int dimension = geometrySpec.getDimension();
    if (dimension == 1) {
        throw new SolverException(simTask.getSimulation().getSolverTaskDescription().getSolverDescription().getDisplayLabel() + " is only supported for simulations with 2D or 3D geometry.");
    }
    Simulation simulation = getSimulationTask().getSimulation();
    SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
    ChomboSolverSpec chomboSolverSpec = solverTaskDescription.getChomboSolverSpec();
    printWriter.println(FVInputFileKeyword.CHOMBO_SPEC_BEGIN);
    printWriter.println(FVInputFileKeyword.DIMENSION + " " + geometrySpec.getDimension());
    Extent extent = geometrySpec.getExtent();
    Origin origin = geometrySpec.getOrigin();
    ISize isize = simulation.getMeshSpecification().getSamplingSize();
    switch(geometrySpec.getDimension()) {
        case 2:
            printWriter.println(FVInputFileKeyword.MESH_SIZE + " " + isize.getX() + " " + isize.getY());
            printWriter.println(FVInputFileKeyword.DOMAIN_SIZE + " " + extent.getX() + " " + extent.getY());
            printWriter.println(FVInputFileKeyword.DOMAIN_ORIGIN + " " + origin.getX() + " " + origin.getY());
            break;
        case 3:
            printWriter.println(FVInputFileKeyword.MESH_SIZE + " " + isize.getX() + " " + isize.getY() + " " + isize.getZ());
            printWriter.println(FVInputFileKeyword.DOMAIN_SIZE + " " + extent.getX() + " " + extent.getY() + " " + extent.getZ());
            printWriter.println(FVInputFileKeyword.DOMAIN_ORIGIN + " " + origin.getX() + " " + origin.getY() + " " + origin.getZ());
            break;
    }
    List<CompartmentSubDomain> featureList = new ArrayList<CompartmentSubDomain>();
    Enumeration<SubDomain> enum1 = simulation.getMathDescription().getSubDomains();
    while (enum1.hasMoreElements()) {
        SubDomain sd = enum1.nextElement();
        if (sd instanceof CompartmentSubDomain) {
            featureList.add((CompartmentSubDomain) sd);
        }
    }
    int numFeatures = featureList.size();
    CompartmentSubDomain[] features = featureList.toArray(new CompartmentSubDomain[0]);
    int[] phases = new int[numFeatures];
    Arrays.fill(phases, -1);
    phases[numFeatures - 1] = 0;
    int[] numAssigned = new int[] { 1 };
    assignPhases(features, numFeatures - 1, phases, numAssigned);
    Map<String, Integer> subDomainPhaseMap = new HashMap<String, Integer>();
    for (int i = 0; i < phases.length; ++i) {
        if (phases[i] == -1) {
            throw new SolverException("Failed to assign a phase to CompartmentSubdomain '" + features[i].getName() + "'. It might be caused by too coarsh a mesh.");
        }
        subDomainPhaseMap.put(features[i].getName(), phases[i]);
    }
    SubVolume[] subVolumes = geometrySpec.getSubVolumes();
    if (geometrySpec.hasImage()) {
        Geometry geometry = (Geometry) BeanUtils.cloneSerializable(simulation.getMathDescription().getGeometry());
        Geometry simGeometry = geometry;
        VCImage img = geometry.getGeometrySpec().getImage();
        int factor = Math.max(Math.max(img.getNumX(), img.getNumY()), img.getNumZ()) < 512 ? 2 : 1;
        ISize distanceMapMeshSize = new ISize(img.getNumX() * factor, img.getNumY() * factor, img.getNumZ() * factor);
        Vect3d deltaX = null;
        boolean bCellCentered = false;
        double dx = 0.5;
        double dy = 0.5;
        double dz = 0.5;
        int Nx = distanceMapMeshSize.getX();
        int Ny = distanceMapMeshSize.getY();
        int Nz = distanceMapMeshSize.getZ();
        if (dimension == 2) {
            // pad the 2D image with itself in order to obtain a 3D image used to compute the distance map
            // because the distance map algorithm is 3D only (using distance to triangles)
            byte[] oldPixels = img.getPixels();
            byte[] newPixels = new byte[oldPixels.length * 3];
            System.arraycopy(oldPixels, 0, newPixels, 0, oldPixels.length);
            System.arraycopy(oldPixels, 0, newPixels, oldPixels.length, oldPixels.length);
            System.arraycopy(oldPixels, 0, newPixels, oldPixels.length * 2, oldPixels.length);
            double distX = geometry.getExtent().getX() / img.getNumX();
            double distY = geometry.getExtent().getY() / img.getNumY();
            // we set the distance on the z axis to something that makes sense
            double distZ = Math.max(distX, distY);
            Extent newExtent = new Extent(geometry.getExtent().getX(), geometry.getExtent().getY(), distZ * 3);
            VCImage newImage = new VCImageUncompressed(null, newPixels, newExtent, img.getNumX(), img.getNumY(), 3);
            // copy the pixel classes too
            ArrayList<VCPixelClass> newPixelClasses = new ArrayList<VCPixelClass>();
            for (VCPixelClass origPixelClass : geometry.getGeometrySpec().getImage().getPixelClasses()) {
                SubVolume origSubvolume = geometry.getGeometrySpec().getImageSubVolumeFromPixelValue(origPixelClass.getPixel());
                newPixelClasses.add(new VCPixelClass(null, origSubvolume.getName(), origPixelClass.getPixel()));
            }
            newImage.setPixelClasses(newPixelClasses.toArray(new VCPixelClass[newPixelClasses.size()]));
            simGeometry = new Geometry(geometry, newImage);
            Nz = 3;
        }
        GeometrySpec simGeometrySpec = simGeometry.getGeometrySpec();
        Extent simExtent = simGeometrySpec.getExtent();
        dx = simExtent.getX() / (Nx - 1);
        dy = simExtent.getY() / (Ny - 1);
        dz = simExtent.getZ() / (Nz - 1);
        if (Math.abs(dx - dy) > 0.1 * Math.max(dx, dy)) {
            dx = Math.min(dx, dy);
            dy = dx;
            Nx = (int) (simExtent.getX() / dx + 1);
            Ny = (int) (simExtent.getY() / dx + 1);
            if (dimension == 3) {
                dz = dx;
                Nz = (int) (simExtent.getZ() / dx + 1);
            }
        }
        deltaX = new Vect3d(dx, dy, dz);
        // one more point in each direction
        distanceMapMeshSize = new ISize(Nx + 1, Ny + 1, Nz + 1);
        Extent distanceMapExtent = new Extent(simExtent.getX() + dx, simExtent.getY() + dy, simExtent.getZ() + dz);
        simGeometrySpec.setExtent(distanceMapExtent);
        GeometrySurfaceDescription geoSurfaceDesc = simGeometry.getGeometrySurfaceDescription();
        geoSurfaceDesc.setVolumeSampleSize(distanceMapMeshSize);
        geoSurfaceDesc.updateAll();
        VCImage vcImage = RayCaster.sampleGeometry(simGeometry, distanceMapMeshSize, bCellCentered);
        SubvolumeSignedDistanceMap[] distanceMaps = DistanceMapGenerator.computeDistanceMaps(simGeometry, vcImage, bCellCentered);
        if (dimension == 2) {
            distanceMaps = DistanceMapGenerator.extractMiddleSlice(distanceMaps);
        }
        printWriter.println(FVInputFileKeyword.SUBDOMAINS + " " + simGeometrySpec.getNumSubVolumes() + " " + FVInputFileKeyword.DISTANCE_MAP);
        for (int i = 0; i < subVolumes.length; i++) {
            File distanceMapFile = new File(workingDirectory, getSimulationTask().getSimulationJobID() + "_" + subVolumes[i].getName() + DISTANCE_MAP_FILE_EXTENSION);
            writeDistanceMapFile(deltaX, distanceMaps[i], distanceMapFile);
            int phase = subDomainPhaseMap.get(subVolumes[i].getName());
            printWriter.println(subVolumes[i].getName() + " " + phase + " " + distanceMapFile.getAbsolutePath());
        }
    } else {
        printWriter.println(FVInputFileKeyword.SUBDOMAINS + " " + geometrySpec.getNumSubVolumes());
        Expression[] rvachevExps = convertAnalyticGeometryToRvachevFunction(geometrySpec);
        for (int i = 0; i < subVolumes.length; i++) {
            if (subVolumes[i] instanceof AnalyticSubVolume) {
                String name = subVolumes[i].getName();
                int phase = subDomainPhaseMap.get(name);
                printWriter.println(name + " " + phase + " ");
                printWriter.println(FVInputFileKeyword.IF + " " + rvachevExps[i].infix() + ";");
                printWriter.println(FVInputFileKeyword.USER + " " + ((AnalyticSubVolume) subVolumes[i]).getExpression().infix() + ";");
            }
        }
    }
    printWriter.println(FVInputFileKeyword.MAX_BOX_SIZE + " " + chomboSolverSpec.getMaxBoxSize());
    printWriter.println(FVInputFileKeyword.FILL_RATIO + " " + chomboSolverSpec.getFillRatio());
    printWriter.println(FVInputFileKeyword.RELATIVE_TOLERANCE + " " + simulation.getSolverTaskDescription().getErrorTolerance().getRelativeErrorTolerance());
    printWriter.println(FVInputFileKeyword.SAVE_VCELL_OUTPUT + " " + chomboSolverSpec.isSaveVCellOutput());
    printWriter.println(FVInputFileKeyword.SAVE_CHOMBO_OUTPUT + " " + chomboSolverSpec.isSaveChomboOutput());
    printWriter.println(FVInputFileKeyword.ACTIVATE_FEATURE_UNDER_DEVELOPMENT + " " + chomboSolverSpec.isActivateFeatureUnderDevelopment());
    printWriter.println(FVInputFileKeyword.SMALL_VOLFRAC_THRESHOLD + " " + chomboSolverSpec.getSmallVolfracThreshold());
    printWriter.println(FVInputFileKeyword.BLOCK_FACTOR + " " + chomboSolverSpec.getBlockFactor());
    printWriter.println(FVInputFileKeyword.TAGS_GROW + " " + chomboSolverSpec.getTagsGrow());
    // Refinement
    int numLevels = chomboSolverSpec.getNumRefinementLevels();
    // Refinements #Levels ratio 1, ratio 2, etc
    printWriter.print(FVInputFileKeyword.REFINEMENTS + " " + (numLevels + 1));
    List<Integer> ratios = chomboSolverSpec.getRefineRatioList();
    for (int i : ratios) {
        printWriter.print(" " + i);
    }
    // write last refinement ratio, fake
    printWriter.println(" 2");
    // membrane rois
    List<RefinementRoi> memRios = chomboSolverSpec.getMembraneRefinementRois();
    printWriter.println(FVInputFileKeyword.REFINEMENT_ROIS + " " + RoiType.Membrane + " " + memRios.size());
    for (RefinementRoi roi : memRios) {
        if (roi.getRoiExpression() == null) {
            throw new SolverException("ROI expression cannot be null");
        }
        // level tagsGrow ROIexpression
        printWriter.println(roi.getLevel() + " " + roi.getRoiExpression().infix() + ";");
    }
    List<RefinementRoi> volRios = chomboSolverSpec.getVolumeRefinementRois();
    printWriter.println(FVInputFileKeyword.REFINEMENT_ROIS + " " + RoiType.Volume + " " + volRios.size());
    for (RefinementRoi roi : volRios) {
        if (roi.getRoiExpression() == null) {
            throw new SolverException("ROI expression cannot be null");
        }
        printWriter.println(roi.getLevel() + " " + roi.getRoiExpression().infix() + ";");
    }
    printWriter.println(FVInputFileKeyword.VIEW_LEVEL + " " + chomboSolverSpec.getViewLevel());
    printWriter.println(FVInputFileKeyword.CHOMBO_SPEC_END);
    printWriter.println();
}