Examples with RefinementRoi org.vcell.chombo.RefinementRoi used on opensource projects

Example 1 with RefinementRoi

use of org.vcell.chombo.RefinementRoi 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 2 with RefinementRoi

use of org.vcell.chombo.RefinementRoi in project vcell by virtualcell.

the class Xmlproducer method getXML.

private Element getXML(ChomboSolverSpec chomboSolverSpec) {
    Element chomboElement = new Element(XMLTags.ChomboSolverSpec);
    Element maxBoxSize = new Element(XMLTags.MaxBoxSizeTag);
    maxBoxSize.setText(chomboSolverSpec.getMaxBoxSize() + "");
    chomboElement.addContent(maxBoxSize);
    if (!chomboSolverSpec.isFinestViewLevel()) {
        Element viewLevel = new Element(XMLTags.ViewLevelTag);
        viewLevel.setText(chomboSolverSpec.getViewLevel() + "");
        chomboElement.addContent(viewLevel);
    }
    Element fillRatio = new Element(XMLTags.FillRatioTag);
    fillRatio.setText(chomboSolverSpec.getFillRatio() + "");
    chomboElement.addContent(fillRatio);
    List<Integer> ratios = chomboSolverSpec.getRefineRatioList();
    if (ratios.size() > 0) {
        Element ratiosElement = new Element(XMLTags.RefineRatios);
        StringBuilder sb = new StringBuilder();
        for (Integer i : ratios) {
            sb.append("," + i);
        }
        ratiosElement.setText(sb.substring(1) + "");
        chomboElement.addContent(ratiosElement);
    }
    Element saveVCellOutput = new Element(XMLTags.SaveVCellOutput);
    saveVCellOutput.setText(chomboSolverSpec.isSaveVCellOutput() + "");
    chomboElement.addContent(saveVCellOutput);
    Element saveChomboOutput = new Element(XMLTags.SaveChomboOutput);
    saveChomboOutput.setText(chomboSolverSpec.isSaveChomboOutput() + "");
    chomboElement.addContent(saveChomboOutput);
    Element activateFeatureUnderDevelopment = new Element(XMLTags.ActivateFeatureUnderDevelopment);
    activateFeatureUnderDevelopment.setText(chomboSolverSpec.isActivateFeatureUnderDevelopment() + "");
    chomboElement.addContent(activateFeatureUnderDevelopment);
    Element smallVolfracThreshold = new Element(XMLTags.SmallVolfracThreshold);
    smallVolfracThreshold.setText(chomboSolverSpec.getSmallVolfracThreshold() + "");
    chomboElement.addContent(smallVolfracThreshold);
    Element blockFactor = new Element(XMLTags.BlockFactorTag);
    blockFactor.setText(chomboSolverSpec.getBlockFactor() + "");
    chomboElement.addContent(blockFactor);
    Element elementTagsGrow = new Element(XMLTags.TagsGrowTag);
    elementTagsGrow.setText(chomboSolverSpec.getTagsGrow() + "");
    chomboElement.addContent(elementTagsGrow);
    Element timeBoundsTag = new Element(XMLTags.TimeBoundTag);
    for (TimeInterval ti : chomboSolverSpec.getTimeIntervalList()) {
        Element e = new Element(XMLTags.TimeIntervalTag);
        e.setAttribute(XMLTags.StartTimeAttrTag, String.valueOf(ti.getStartingTime()));
        e.setAttribute(XMLTags.EndTimeAttrTag, String.valueOf(ti.getEndingTime()));
        e.setAttribute(XMLTags.TimeStepAttrTag, String.valueOf(ti.getTimeStep()));
        e.setAttribute(XMLTags.OutputTimeStepAttrTag, String.valueOf(ti.getOutputTimeStep()));
        timeBoundsTag.addContent(e);
    }
    chomboElement.addContent(timeBoundsTag);
    Element meshRefinement = new Element(XMLTags.MeshRefinementTag);
    for (RefinementRoi roi : chomboSolverSpec.getRefinementRois()) {
        Element roiElement = new Element(XMLTags.RefinementRoiTag);
        roiElement.setAttribute(XMLTags.RefineRoiLevelAttrTag, String.valueOf(roi.getLevel()));
        roiElement.setAttribute(XMLTags.RefinementRoiTypeAttrTag, roi.getType().name());
        if (roi.getRoiExpression() != null) {
            Element expElement = new Element(XMLTags.ROIExpressionTag);
            expElement.setText(roi.getRoiExpression().infix() + ";");
            roiElement.addContent(expElement);
        }
        meshRefinement.addContent(roiElement);
    }
    chomboElement.addContent(meshRefinement);
    return chomboElement;
}

Example 3 with RefinementRoi

use of org.vcell.chombo.RefinementRoi in project vcell by virtualcell.

the class XmlReader method getChomboSolverSpec.

private ChomboSolverSpec getChomboSolverSpec(SolverTaskDescription solverTaskDesc, Element element, int dimension) throws XmlParseException {
    int maxBoxSize = parseIntWithDefault(element, XMLTags.MaxBoxSizeTag, ChomboSolverSpec.getDefaultMaxBoxSize(dimension));
    double fillRatio = parseDoubleWithDefault(element, XMLTags.FillRatioTag, ChomboSolverSpec.getDefaultFillRatio());
    boolean bSaveVCellOutput = parseBooleanWithDefault(element, XMLTags.SaveVCellOutput, true);
    boolean bSaveChomboOutput = parseBooleanWithDefault(element, XMLTags.SaveChomboOutput, false);
    Element childElement = element.getChild(XMLTags.RefineRatios, vcNamespace);
    List<Integer> refineRatioList = null;
    if (childElement != null) {
        String text = childElement.getText();
        if (text != null && !text.isEmpty()) {
            StringTokenizer st = new StringTokenizer(text, ",");
            if (st.hasMoreTokens()) {
                refineRatioList = new ArrayList<Integer>();
                while (st.hasMoreElements()) {
                    String token = st.nextToken();
                    if (token != null) {
                        int n = Integer.parseInt(token);
                        refineRatioList.add(n);
                    }
                }
            }
        }
    }
    Integer viewLevel = null;
    try {
        ChomboSolverSpec css = new ChomboSolverSpec(maxBoxSize, fillRatio, viewLevel, bSaveVCellOutput, bSaveChomboOutput, refineRatioList);
        double smallVolfracThreshold = parseDoubleWithDefault(element, XMLTags.SmallVolfracThreshold, 0);
        int blockFactor = parseIntWithDefault(element, XMLTags.BlockFactorTag, ChomboSolverSpec.DEFAULT_BLOCK_FACTOR);
        boolean bActivateFeatureUnderDevelopment = parseBooleanWithDefault(element, XMLTags.ActivateFeatureUnderDevelopment, false);
        css.setSmallVolfracThreshold(smallVolfracThreshold);
        css.setActivateFeatureUnderDevelopment(bActivateFeatureUnderDevelopment);
        css.setBlockFactor(blockFactor);
        int tagsGrow = parseIntWithDefault(element, XMLTags.TagsGrowTag, ChomboSolverSpec.defaultTagsGrow);
        css.setTagsGrow(tagsGrow);
        Element timeBoundsElement = element.getChild(XMLTags.TimeBoundTag, vcNamespace);
        List<Element> timeIntervalElementList = null;
        boolean noTimeBounds = false;
        if (timeBoundsElement == null) {
            noTimeBounds = true;
        } else {
            timeIntervalElementList = timeBoundsElement.getChildren(XMLTags.TimeIntervalTag, vcNamespace);
            if (timeIntervalElementList.size() == 0) {
                noTimeBounds = true;
            }
        }
        if (noTimeBounds) {
            // old format
            double startTime = 0;
            double endTime = solverTaskDesc.getTimeBounds().getEndingTime();
            double timeStep = solverTaskDesc.getTimeStep().getDefaultTimeStep();
            double outputTimeStep = ((UniformOutputTimeSpec) solverTaskDesc.getOutputTimeSpec()).getOutputTimeStep();
            try {
                TimeInterval ti = new TimeInterval(startTime, endTime, timeStep, outputTimeStep);
                css.addTimeInterval(ti);
            } catch (IllegalArgumentException ex) {
                css.addTimeInterval(TimeInterval.getDefaultTimeInterval());
            }
        } else {
            for (Element e : timeIntervalElementList) {
                String s = e.getAttributeValue(XMLTags.StartTimeAttrTag);
                double startTime = Double.valueOf(s);
                s = e.getAttributeValue(XMLTags.EndTimeAttrTag);
                double endTime = Double.valueOf(s);
                s = e.getAttributeValue(XMLTags.TimeStepAttrTag);
                double timeStep = Double.valueOf(s);
                s = e.getAttributeValue(XMLTags.OutputTimeStepAttrTag);
                double outputTimeStep = Double.valueOf(s);
                TimeInterval ti = new TimeInterval(startTime, endTime, timeStep, outputTimeStep);
                css.addTimeInterval(ti);
            }
        }
        Element meshRefineElement = element.getChild(XMLTags.MeshRefinementTag, vcNamespace);
        if (meshRefineElement != null) {
            if (meshRefineElement.getChildren().size() != 0) {
                // in old model, if there is no refinement, set view level to finest
                // only set viewLevel when meshRefinement has children
                Element viewLevelChild = element.getChild(XMLTags.ViewLevelTag, vcNamespace);
                if (viewLevelChild != null) {
                    viewLevel = parseIntWithDefault(element, XMLTags.ViewLevelTag, 0);
                    css.setViewLevel(viewLevel);
                }
            }
            List<Element> levelElementList = meshRefineElement.getChildren(XMLTags.RefinementRoiTag, vcNamespace);
            for (Element levelElement : levelElementList) {
                String levelStr = levelElement.getAttributeValue(XMLTags.RefineRoiLevelAttrTag);
                int level = 1;
                if (levelStr != null) {
                    level = Integer.parseInt(levelStr);
                }
                String type = levelElement.getAttributeValue(XMLTags.RefinementRoiTypeAttrTag);
                RoiType roiType = RoiType.Membrane;
                if (type != null) {
                    try {
                        roiType = RoiType.valueOf(type);
                    } catch (Exception ex) {
                    // ignore
                    }
                }
                Element expElement = levelElement.getChild(XMLTags.ROIExpressionTag, vcNamespace);
                String roiExp = null;
                if (expElement != null) {
                    roiExp = expElement.getText();
                    RefinementRoi roi = new RefinementRoi(roiType, level, roiExp);
                    css.addRefinementRoi(roi);
                }
            }
        }
        return css;
    } catch (ExpressionException e) {
        throw new XmlParseException(e);
    }
}