Examples with Vect3d cbit.vcell.render.Vect3d used on opensource projects

Example 6 with Vect3d

use of cbit.vcell.render.Vect3d in project vcell by virtualcell.

the class CSGObjectPropertiesPanel method stringToVect3d.

static Vect3d stringToVect3d(String str) {
    StringTokenizer st = new StringTokenizer(str, ", ");
    double x = Double.parseDouble(st.nextToken());
    double y = Double.parseDouble(st.nextToken());
    double z = Double.parseDouble(st.nextToken());
    Vect3d vect3d = new Vect3d(x, y, z);
    return vect3d;
}

Example 7 with Vect3d

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

Example 8 with Vect3d

use of cbit.vcell.render.Vect3d in project vcell by virtualcell.

the class CSGScale method updateTransform.

private void updateTransform() {
    Affine forward = new Affine();
    forward.setScale(scale);
    Affine inverse = new Affine();
    inverse.setScale(new Vect3d(1.0 / scale.getX(), 1.0 / scale.getY(), 1.0 / scale.getZ()));
    setTransforms(forward, inverse);
}

Example 9 with Vect3d

use of cbit.vcell.render.Vect3d 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 10 with Vect3d

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