Examples with Triangle cbit.vcell.geometry.surface.Triangle used on opensource projects

Example 1 with Triangle

use of cbit.vcell.geometry.surface.Triangle 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 Triangle

use of cbit.vcell.geometry.surface.Triangle 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 Triangle

use of cbit.vcell.geometry.surface.Triangle in project vcell by virtualcell.

the class SphereTestFixture method getTestSphere.

/**
 * @return test spherical tessellation
 * @throws ProgrammingException if not {@link #isSphereTest()}
 */
public List<? extends Polygon> getTestSphere() {
    if (isSphereTest()) {
        IcoSphere icoSphere = IcoSphere.get();
        List<PolygonImmutable> tessel = icoSphere.getTessellation(level);
        ThreeSpacePoint centerPoint = new ThreeSpacePointImmutable(center[0], center[1], center[2]);
        List<Triangle> tl = Triangle.scale(tessel, radius, centerPoint);
        return tl;
    }
    throw new ProgrammingException("getTestSphere called on non test Simulation");
}

Example 4 with Triangle

use of cbit.vcell.geometry.surface.Triangle in project vcell by virtualcell.

the class SmoldynFileWriter method writeInitialConcentration.

private int writeInitialConcentration(ParticleInitialConditionConcentration initialConcentration, SubDomain subDomain, Variable variable, String variableName, StringBuilder sb) throws ExpressionException, MathException {
    SimpleSymbolTable simpleSymbolTable = new SimpleSymbolTable(new String[] { ReservedVariable.X.getName(), ReservedVariable.Y.getName(), ReservedVariable.Z.getName() });
    Expression disExpression = new Expression(initialConcentration.getDistribution());
    disExpression.bindExpression(simulationSymbolTable);
    disExpression = simulationSymbolTable.substituteFunctions(disExpression).flatten();
    disExpression.bindExpression(simpleSymbolTable);
    double[] values = new double[3];
    if (dimension == 1) {
        if (disExpression.getSymbolBinding(ReservedVariable.Y.getName()) != null || disExpression.getSymbolBinding(ReservedVariable.Z.getName()) != null) {
            throw new MathException(VCellErrorMessages.getSmoldynWrongCoordinates("'y' or 'z'", dimension, variable, disExpression));
        }
    } else if (dimension == 2) {
        if (disExpression.getSymbolBinding(ReservedVariable.Z.getName()) != null) {
            throw new MathException(VCellErrorMessages.getSmoldynWrongCoordinates("'z'", dimension, variable, disExpression));
        }
    }
    int totalCount = 0;
    StringBuilder localsb = new StringBuilder();
    if (subDomain instanceof CompartmentSubDomain) {
        MeshSpecification meshSpecification = simulation.getMeshSpecification();
        ISize sampleSize = meshSpecification.getSamplingSize();
        int numX = sampleSize.getX();
        int numY = dimension < 2 ? 1 : sampleSize.getY();
        int numZ = dimension < 3 ? 1 : sampleSize.getZ();
        boolean bCellCentered = simulation.hasCellCenteredMesh();
        double dx = meshSpecification.getDx(bCellCentered);
        double dy = meshSpecification.getDy(bCellCentered);
        double dz = meshSpecification.getDz(bCellCentered);
        Origin origin = resampledGeometry.getGeometrySpec().getOrigin();
        double ox = origin.getX();
        double oy = origin.getY();
        double oz = origin.getZ();
        Extent extent = resampledGeometry.getExtent();
        double ex = extent.getX();
        double ey = extent.getY();
        double ez = extent.getZ();
        int offset = 0;
        for (int k = 0; k < numZ; k++) {
            double centerz = oz + k * dz;
            double loz = Math.max(oz, centerz - dz / 2);
            double hiz = Math.min(oz + ez, centerz + dz / 2);
            double lz = hiz - loz;
            values[2] = centerz;
            for (int j = 0; j < numY; j++) {
                double centery = oy + j * dy;
                double loy = Math.max(oy, centery - dy / 2);
                double hiy = Math.min(oy + ey, centery + dy / 2);
                values[1] = centery;
                double ly = hiy - loy;
                for (int i = 0; i < numX; i++) {
                    int regionIndex = resampledGeometry.getGeometrySurfaceDescription().getRegionImage().getRegionInfoFromOffset(offset).getRegionIndex();
                    offset++;
                    GeometricRegion region = resampledGeometry.getGeometrySurfaceDescription().getGeometricRegions(regionIndex);
                    if (region instanceof VolumeGeometricRegion) {
                        if (!((VolumeGeometricRegion) region).getSubVolume().getName().equals(subDomain.getName())) {
                            continue;
                        }
                    }
                    double centerx = ox + i * dx;
                    double lox = Math.max(ox, centerx - dx / 2);
                    double hix = Math.min(ox + ex, centerx + dx / 2);
                    double lx = hix - lox;
                    values[0] = centerx;
                    double volume = lx;
                    if (dimension > 1) {
                        volume *= ly;
                        if (dimension > 2) {
                            volume *= lz;
                        }
                    }
                    double expectedCount = disExpression.evaluateVector(values) * volume;
                    if (expectedCount <= 0) {
                        continue;
                    }
                    long count = dist.nextPoisson(expectedCount);
                    if (count <= 0) {
                        continue;
                    }
                    totalCount += count;
                    localsb.append(SmoldynVCellMapper.SmoldynKeyword.mol + " " + count + " " + variableName + " " + (float) lox + "-" + (float) hix);
                    if (lg.isDebugEnabled()) {
                        lg.debug("Component subdomain " + variableName + " count " + count);
                    }
                    if (dimension > 1) {
                        localsb.append(" " + loy + "-" + hiy);
                        if (dimension > 2) {
                            localsb.append(" " + loz + "-" + hiz);
                        }
                    }
                    localsb.append("\n");
                }
            }
        }
        // otherwise we append the distributed molecules in different small boxes
        try {
            subsituteFlattenToConstant(disExpression);
            sb.append(SmoldynVCellMapper.SmoldynKeyword.compartment_mol);
            sb.append(" " + totalCount + " " + variableName + " " + subDomain.getName() + "\n");
        } catch (// can not be evaluated to a constant
        Exception e) {
            sb.append(localsb);
        }
    } else if (subDomain instanceof MembraneSubDomain) {
        ArrayList<TrianglePanel> trianglePanelList = membraneSubdomainTriangleMap.get(subDomain);
        for (TrianglePanel trianglePanel : trianglePanelList) {
            Triangle triangle = trianglePanel.triangle;
            switch(dimension) {
                case 1:
                    values[0] = triangle.getNodes(0).getX();
                    break;
                case 2:
                    {
                        double centroidX = triangle.getNodes(0).getX();
                        double centroidY = triangle.getNodes(0).getY();
                        if (triangle.getNodes(0).getX() == triangle.getNodes(1).getX() && triangle.getNodes(0).getY() == triangle.getNodes(1).getY()) {
                            centroidX += triangle.getNodes(2).getX();
                            centroidY += triangle.getNodes(2).getY();
                        } else {
                            centroidX += triangle.getNodes(1).getX();
                            centroidY += triangle.getNodes(1).getY();
                        }
                        values[0] = centroidX / 2;
                        values[1] = centroidY / 2;
                        break;
                    }
                case 3:
                    {
                        double centroidX = triangle.getNodes(0).getX() + triangle.getNodes(1).getX() + triangle.getNodes(2).getX();
                        double centroidY = triangle.getNodes(0).getY() + triangle.getNodes(1).getY() + triangle.getNodes(2).getY();
                        double centroidZ = triangle.getNodes(0).getZ() + triangle.getNodes(1).getZ() + triangle.getNodes(2).getZ();
                        values[0] = centroidX / 3;
                        values[1] = centroidY / 3;
                        values[2] = centroidZ / 3;
                        break;
                    }
            }
            double expectedCount = disExpression.evaluateVector(values) * triangle.getArea();
            if (expectedCount <= 0) {
                continue;
            }
            long count = dist.nextPoisson(expectedCount);
            if (count <= 0) {
                continue;
            }
            totalCount += count;
            if (lg.isDebugEnabled()) {
                lg.debug("Membrane subdomain " + subDomain.getName() + ' ' + variableName + " count " + count);
            }
            localsb.append(SmoldynVCellMapper.SmoldynKeyword.surface_mol + " " + count + " " + variableName + " " + subDomain.getName() + " " + SmoldynVCellMapper.SmoldynKeyword.tri + " " + trianglePanel.name + "\n");
        }
        // otherwise we append the distributed molecules in different small boxes
        try {
            subsituteFlattenToConstant(disExpression);
            sb.append(SmoldynVCellMapper.SmoldynKeyword.surface_mol);
            sb.append(" " + totalCount + " " + variableName + " " + subDomain.getName() + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.all + "\n");
        } catch (// can not be evaluated to a constant
        Exception e) {
            sb.append(localsb);
        }
    }
    if (lg.isDebugEnabled()) {
        lg.debug("Subdomain " + subDomain.getName() + ' ' + variableName + " total count " + totalCount);
    }
    return totalCount;
}