Examples with Coordinate org.vcell.util.Coordinate used on opensource projects

Example 76 with Coordinate

use of org.vcell.util.Coordinate in project vcell by virtualcell.

the class SpatialSelectionVolume method sampleSymmetric.

/**
 * Insert the method's description here.
 * Creation date: (10/5/2004 7:19:49 AM)
 */
private SSHelper sampleSymmetric() throws Exception {
    // 
    if (!isSymmetric()) {
        return null;
    }
    SampledCurve ssCurve = getCurveSelectionInfo().getCurve().getSampledCurve();
    Vector pointsV = ssCurve.getControlPointsVector();
    if (pointsV.size() < 2) {
        return null;
    }
    final int INDEX_SPACE_INCR = 100;
    int[] indexes = new int[INDEX_SPACE_INCR];
    int indexCounter = 0;
    Vector snappedV = new Vector();
    CoordinateIndex lastCI = null;
    for (int i = 0; i < pointsV.size(); i += 1) {
        if (i > 0) {
            CoordinateIndex c1 = getMesh().getCoordinateIndexFromFractionalIndex(getMesh().getFractionalCoordinateIndex((Coordinate) pointsV.elementAt(i - 1)));
            CoordinateIndex c2 = getMesh().getCoordinateIndexFromFractionalIndex(getMesh().getFractionalCoordinateIndex((Coordinate) pointsV.elementAt(i)));
            int dx = c2.x - c1.x;
            if (dx != 0) {
                dx /= Math.abs(dx);
            }
            int dy = c2.y - c1.y;
            if (dy != 0) {
                dy /= Math.abs(dy);
            }
            int dz = c2.z - c1.z;
            if (dz != 0) {
                dz /= Math.abs(dz);
            }
            while (true) {
                if (!c1.compareEqual(lastCI)) {
                    if (indexCounter == indexes.length) {
                        // Make more space
                        int[] temp = new int[indexes.length + INDEX_SPACE_INCR];
                        System.arraycopy(indexes, 0, temp, 0, indexes.length);
                        indexes = temp;
                    }
                    indexes[indexCounter] = getConvertedIndexFromCI(c1);
                    snappedV.add(getMesh().getCoordinate(c1));
                    indexCounter += 1;
                }
                if (c1.compareEqual(c2)) {
                    break;
                }
                c1.x += dx;
                c1.y += dy;
                c1.z += dz;
                // sanity check to prevent infinite loop, this shouldn't happen
                if (c1.x < 0 || c1.y < 0 || c1.z < 0 || c1.x >= getMesh().getSizeX() || c1.y >= getMesh().getSizeY() || c1.z >= getMesh().getSizeZ()) {
                    throw new Exception(this.getClass().getName() + ".sampleSymmetric failed");
                }
            }
            ;
            lastCI = c2;
        }
    }
    if (indexCounter > 0) {
        return makeSSHelper(indexes, indexCounter, snappedV, false);
    }
    return null;
}

Example 77 with Coordinate

use of org.vcell.util.Coordinate in project vcell by virtualcell.

the class SpatialSelectionVolume method isSymmetric.

/**
 * Insert the method's description here.
 * Creation date: (10/3/2004 12:54:59 PM)
 * @return boolean
 */
private boolean isSymmetric() {
    if (!(getCurveSelectionInfo().getCurve() instanceof SampledCurve)) {
        return false;
    }
    SampledCurve ssCurve = getCurveSelectionInfo().getCurve().getSampledCurve();
    Vector pointsV = ssCurve.getControlPointsVector();
    if (pointsV.size() == 1) {
        return true;
    }
    for (int i = 0; i < pointsV.size(); i += 1) {
        if (i > 0) {
            CoordinateIndex c1 = getMesh().getCoordinateIndexFromFractionalIndex(getMesh().getFractionalCoordinateIndex((Coordinate) pointsV.elementAt(i - 1)));
            CoordinateIndex c2 = getMesh().getCoordinateIndexFromFractionalIndex(getMesh().getFractionalCoordinateIndex((Coordinate) pointsV.elementAt(i)));
            int dx = Math.abs(c1.x - c2.x);
            int dy = Math.abs(c1.y - c2.y);
            int dz = Math.abs(c1.z - c2.z);
            int symCount = 0;
            if (dx != 0 && dy != 0 && dx != dy) {
                symCount += 1;
            }
            if (dx != 0 && dz != 0 && dx != dz) {
                symCount += 1;
            }
            if (dy != 0 && dz != 0 && dy != dz) {
                symCount += 1;
            }
            if (symCount != 0) {
                return false;
            }
        }
    }
    return true;
}

Example 78 with Coordinate

use of org.vcell.util.Coordinate in project vcell by virtualcell.

the class SpatialSelectionVolume method getIndexSamples.

/**
 * Insert the method's description here.
 * Creation date: (2/26/2001 6:17:10 PM)
 * @return int[]
 * @param numSamples int
 */
public SSHelper getIndexSamples(double begin, double end) {
    if (getCurveSelectionInfo().getCurve() instanceof SinglePoint) {
        return new SSHelper(new Coordinate[] { getCurveSelectionInfo().getCurve().getCoordinate(0) }, new int[] { getConvertedIndexFromWC(getCurveSelectionInfo().getCurve().getCoordinate(0)) }, getVariableType(), new int[] { -1 });
    }
    // Try simple sampling
    SSHelper ssvHelper = null;
    if (begin == 0.0 && end == 1.0) {
        try {
            ssvHelper = sampleSymmetric();
        } catch (Throwable e) {
        // Do nothing, we'll try sampling below
        }
    }
    if (ssvHelper != null) {
        return ssvHelper;
    }
    // 
    // Find continuous mesh elements along the curve
    // 
    final int SAMPLE_MULT = 10;
    // Determine reasonable sample size
    double smallestMeshCellDimension = getSmallestMeshCellDimensionLength();
    double uSpatialLength = getCurveSelectionInfo().getCurve().getSpatialLength() * (Math.abs(end - begin));
    int SSV_NUM_SAMPLES = (int) (Math.ceil(uSpatialLength / smallestMeshCellDimension) * (double) SAMPLE_MULT);
    if (SSV_NUM_SAMPLES < 2) {
        SSV_NUM_SAMPLES = 2;
    }
    // 
    double delta = (end - begin) / (double) SSV_NUM_SAMPLES;
    int[] uniquePoints = new int[SSV_NUM_SAMPLES];
    int uniquePointCount = 0;
    Vector<Coordinate> wcV = new Vector<Coordinate>();
    // Vector ciV = new Vector();
    // double[] cellStepArr = new double[SAMPLE_MULT*10];//big enough
    // big enough
    double[] cellStepArr = new double[SSV_NUM_SAMPLES];
    int cellStepCounter = 0;
    int lastISample = -1;
    int currISample;
    double currentU = begin;
    for (int index = 0; index < SSV_NUM_SAMPLES; index += 1) {
        boolean isLastLoop = index == (SSV_NUM_SAMPLES - 1);
        if (isLastLoop) {
            currISample = getMeshIndexFromU(end);
            currentU = end;
        } else {
            currISample = getMeshIndexFromU(currentU);
        }
        if ((index == 0) || (lastISample != currISample)) {
            uniquePoints[uniquePointCount] = getIndex(currentU);
            double midU = midpoint(cellStepArr, cellStepCounter);
            if (index == 0) {
                midU = begin;
            }
            if (isLastLoop) {
                midU = end;
            }
            Coordinate wc = getSamplingWorldCoordinate(midU);
            if (index == 0 || isLastLoop || (uniquePointCount > 1)) {
                wcV.add(wc);
            }
            // CoordinateIndex ci = getCoordinateIndexFromWC(wc);
            // ciV.add(ci);
            uniquePointCount += 1;
            // if(ciV.size() > 1 &&!areTouching(((CoordinateIndex)ciV.elementAt(ciV.size()-1)),((CoordinateIndex)ciV.elementAt(ciV.size()-2)))){
            // //this shouldn't happen if our sampling is fine enough
            // //all sampled mesh elements should be "touching"
            // }
            cellStepCounter = 0;
        } else if (isLastLoop && (wcV.size() != uniquePointCount)) {
            wcV.add(getSamplingWorldCoordinate(end));
        }
        cellStepArr[cellStepCounter] = currentU;
        cellStepCounter += 1;
        currentU += delta;
        lastISample = currISample;
    }
    if (uniquePointCount > 0) {
        return makeSSHelper(uniquePoints, uniquePointCount, wcV, true);
    }
    // This shouldn't happen
    throw new RuntimeException(this.getClass().getName() + " couldn't generate sampling");
}

Example 79 with Coordinate

use of org.vcell.util.Coordinate in project vcell by virtualcell.

the class SpatialSelectionVolume method resampleMeshBoundaries.

/**
 * Insert the method's description here.
 * Creation date: (10/6/2004 11:31:20 AM)
 */
private SSHelper resampleMeshBoundaries(int[] indexes, Coordinate[] wcV, boolean bRecenter) throws Exception {
    if (wcV.length < 2) {
        return null;
    }
    final int MEMBRANE_BOUNDARY = -1;
    final int NON_MEMBRANE_BOUNDARY = -2;
    final int MEMB_INDEX = 0;
    final int IN_VOL = 1;
    final int OUT_VOL = 2;
    int newCount = 0;
    int[] newIndexes = new int[indexes.length * 3];
    int[][][] newCrossingMembraneIndex = new int[newIndexes.length][4][3];
    Coordinate[] newWCV = new Coordinate[wcV.length * 3];
    for (int i = 1; i < wcV.length; i += 1) {
        CoordinateIndex ci1 = getCoordinateIndexFromWC(wcV[i]);
        CoordinateIndex ci2 = getCoordinateIndexFromWC(wcV[i - 1]);
        CoordinateIndex ci3 = null;
        CoordinateIndex ci4 = null;
        CoordinateIndex[][] faces = null;
        if (areTouchingFace(ci1, ci2)) {
            faces = new CoordinateIndex[][] { { ci1, ci2 } };
        } else if (areTouching(ci1, ci2)) {
            // must be a corner
            // determine the 4 faces that coincide with this corner
            int dx = ci2.x - ci1.x;
            int dy = ci2.y - ci1.y;
            int dz = ci2.z - ci1.z;
            if (dx == 0 && dy != 0 && dz != 0) {
                ci3 = new CoordinateIndex(ci1.x, ci1.y + dy, ci1.z);
                ci4 = new CoordinateIndex(ci1.x, ci1.y, ci1.z + dz);
            } else if (dx != 0 && dy == 0 && dz != 0) {
                ci3 = new CoordinateIndex(ci1.x + dx, ci1.y, ci1.z);
                ci4 = new CoordinateIndex(ci1.x, ci1.y, ci1.z + dz);
            } else if (dx != 0 && dy != 0 && dz == 0) {
                ci3 = new CoordinateIndex(ci1.x, ci1.y + dy, ci1.z);
                ci4 = new CoordinateIndex(ci1.x + dx, ci1.y, ci1.z);
            } else {
                throw new Exception(this.getClass().getName() + ".resampleMeshBoundaries Couldn't adjust for corners");
            }
            faces = new CoordinateIndex[][] { { ci1, ci3 }, { ci1, ci4 }, { ci2, ci3 }, { ci2, ci4 } };
        } else {
            throw new Exception(this.getClass().getName() + ".resampleMeshBoundaries Expecting touching elements");
        }
        // System.out.println("p1="+i+" p2="+(i-1));
        boolean[] bCrossMembraneArr = new boolean[faces.length];
        java.util.Arrays.fill(bCrossMembraneArr, false);
        int[][] crossMembraneIndexInOutArr = new int[bCrossMembraneArr.length][3];
        Coordinate[] intersectArr = new Coordinate[faces.length];
        int intersectNotNullCount = 0;
        // Find out where the line segment between volume elements intersect the face(s)
        for (int j = 0; j < faces.length; j += 1) {
            crossMembraneIndexInOutArr[j][MEMB_INDEX] = -1;
            crossMembraneIndexInOutArr[j][IN_VOL] = -1;
            crossMembraneIndexInOutArr[j][OUT_VOL] = -1;
            intersectArr[j] = lineMeshFaceIntersect3D(getMesh().getCoordinate(faces[j][0]), getMesh().getCoordinate(faces[j][1]), wcV[i], wcV[i - 1]);
            intersectNotNullCount += (intersectArr[j] != null ? 1 : 0);
            if (getMesh().getMembraneElements() != null) {
                // Find out if this face is part of a Membrane
                int vi1 = getMesh().getVolumeIndex(faces[j][0]);
                int vi2 = getMesh().getVolumeIndex(faces[j][1]);
                for (int k = 0; k < getMesh().getMembraneElements().length; k += 1) {
                    int inside = getMesh().getMembraneElements()[k].getInsideVolumeIndex();
                    int outside = getMesh().getMembraneElements()[k].getOutsideVolumeIndex();
                    if ((vi1 == inside && vi2 == outside) || (vi1 == outside && vi2 == inside)) {
                        bCrossMembraneArr[j] = true;
                        crossMembraneIndexInOutArr[j][MEMB_INDEX] = k;
                        crossMembraneIndexInOutArr[j][IN_VOL] = inside;
                        crossMembraneIndexInOutArr[j][OUT_VOL] = outside;
                        break;
                    }
                }
            }
        // System.out.println("----- vi1="+vi1+" vi2="+vi2+" face="+j+" intersect="+intersectArr[j]+" bCrossMembrane="+bCrossMembraneArr[j]);
        }
        // we will use them later to determing the midpoint to adjust our data sample WorldCoordinates
        if (i == 1) {
            newIndexes[newCount] = indexes[0];
            newWCV[newCount] = wcV[0];
            newCount += 1;
        }
        if (faces.length == 1 && intersectNotNullCount == 1) {
            // through face, use where it intersects and mark if has Membrane
            newIndexes[newCount] = (bCrossMembraneArr[0] ? MEMBRANE_BOUNDARY : NON_MEMBRANE_BOUNDARY);
            newCrossingMembraneIndex[newCount] = crossMembraneIndexInOutArr;
            newWCV[newCount] = intersectArr[0];
            newCount += 1;
        } else if (faces.length == 4 && intersectNotNullCount == 4) {
            // Through a corner, calculate the corner Coordinate and mark if has Membrane
            Coordinate mc1 = getMesh().getCoordinate(ci1);
            Coordinate mc2 = getMesh().getCoordinate(ci2);
            Coordinate mc3 = getMesh().getCoordinate(ci3);
            Coordinate mc4 = getMesh().getCoordinate(ci4);
            Coordinate corner = new Coordinate((mc1.getX() + mc2.getX() + mc3.getX() + mc4.getX()) / 4.0, (mc1.getY() + mc2.getY() + mc3.getY() + mc4.getY()) / 4.0, (mc1.getZ() + mc2.getZ() + mc3.getZ() + mc4.getZ()) / 4.0);
            newIndexes[newCount] = (bCrossMembraneArr[0] || bCrossMembraneArr[1] || bCrossMembraneArr[2] || bCrossMembraneArr[3] ? MEMBRANE_BOUNDARY : NON_MEMBRANE_BOUNDARY);
            // newCrossingMembraneIndex[newCount] = crossMembraneIndexInOutArr;
            newWCV[newCount] = corner;
            newCount += 1;
        } else {
            throw new Exception(this.getClass().getName() + " Unexpected intersection result");
        }
        newIndexes[newCount] = indexes[i];
        newWCV[newCount] = wcV[i];
        newCount += 1;
    }
    // Count how many points we will finally end up with after markers are removed and
    // including the points that will be added for membrane intersections
    // also calculate the midpoints of the original data between the intersect points
    // first and last
    int finalCount = 2;
    // Adjust distances to center between boundary markers (inersections)
    for (int i = 0; i < newCount; i += 1) {
        if (newIndexes[i] == MEMBRANE_BOUNDARY || newIndexes[i] == NON_MEMBRANE_BOUNDARY) {
            if ((i + 2) < newCount) {
                // if(bRecenter){
                // Coordinate b0 = newWCV[i];
                // Coordinate b1 = newWCV[i+2];
                // newWCV[i+1] = new Coordinate((b0.getX()+b1.getX())/2.0,(b0.getY()+b1.getY())/2.0,(b0.getZ()+b1.getZ())/2.0);
                // }else{
                // TODO: why is this here, is this really what we want? ... A = A;
                newWCV[i + 1] = newWCV[i + 1];
                // }
                finalCount += 1;
            }
        }
        if (newIndexes[i] == MEMBRANE_BOUNDARY) {
            finalCount += 2;
        }
    }
    // Final pass, remove markers, add additional points at membrane intersections
    int[] finalIndexes = new int[finalCount];
    int[] finalCrossingMembrane = null;
    Coordinate[] finalWC = new Coordinate[finalCount];
    finalCount = 0;
    for (int i = 0; i < newCount; i += 1) {
        if (newIndexes[i] != MEMBRANE_BOUNDARY && newIndexes[i] != NON_MEMBRANE_BOUNDARY) {
            // Get all original data points
            finalIndexes[finalCount] = newIndexes[i];
            finalWC[finalCount] = newWCV[i];
            finalCount += 1;
        } else if (newIndexes[i] == MEMBRANE_BOUNDARY) {
            for (int j = 0; j < newCrossingMembraneIndex[i].length; j++) {
                if (newCrossingMembraneIndex[i][j][MEMB_INDEX] != -1 && ((newCrossingMembraneIndex[i][j][IN_VOL] == newIndexes[i - 1] && newCrossingMembraneIndex[i][j][OUT_VOL] == newIndexes[i + 1]) || (newCrossingMembraneIndex[i][j][OUT_VOL] == newIndexes[i - 1] && newCrossingMembraneIndex[i][j][IN_VOL] == newIndexes[i + 1]))) {
                    if (finalCrossingMembrane == null) {
                        // make once only if needed
                        finalCrossingMembrane = new int[finalIndexes.length];
                        Arrays.fill(finalCrossingMembrane, -1);
                    }
                    finalCrossingMembrane[finalCount] = newCrossingMembraneIndex[i][j][MEMB_INDEX];
                    finalCrossingMembrane[finalCount + 1] = newCrossingMembraneIndex[i][j][MEMB_INDEX];
                    break;
                }
            }
            // Add 2 new membrane intersection points, one for each side of the face
            finalIndexes[finalCount] = newIndexes[i - 1];
            // finalWC[finalCount] = offsetCoordinate(newWCV[i],newWCV[i-1]);
            finalWC[finalCount] = newWCV[i];
            finalCount += 1;
            finalIndexes[finalCount] = newIndexes[i + 1];
            // finalWC[finalCount] = offsetCoordinate(newWCV[i],newWCV[i+1]);
            finalWC[finalCount] = newWCV[i];
            finalCount += 1;
        }
    }
    return new SSHelper(finalWC, finalIndexes, getVariableType(), finalCrossingMembrane);
}

Example 80 with Coordinate

use of org.vcell.util.Coordinate in project vcell by virtualcell.

the class SmoldynFileWriter method writeWallSurfaces.

private void writeWallSurfaces() throws SolverException {
    GeometrySurfaceDescription geometrySurfaceDescription = resampledGeometry.getGeometrySurfaceDescription();
    GeometrySpec geometrySpec = resampledGeometry.getGeometrySpec();
    SubVolume[] subVolumes = geometrySpec.getSubVolumes();
    printWriter.println("# boundaries");
    Origin origin = geometrySpec.getOrigin();
    Extent extent = geometrySpec.getExtent();
    Coordinate lowWall = new Coordinate(origin.getX(), origin.getY(), origin.getZ());
    Coordinate highWall = new Coordinate(origin.getX() + extent.getX(), origin.getY() + extent.getY(), origin.getZ() + extent.getZ());
    // potential artifact.
    if (bHasNoSurface) {
        SubDomain subDomain0 = mathDesc.getSubDomains().nextElement();
        CompartmentSubDomain compartSubDomain0 = null;
        compartSubDomain0 = (CompartmentSubDomain) subDomain0;
        // x
        if (compartSubDomain0.getBoundaryConditionXm().isPERIODIC()) {
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 0 " + lowWall.getX() + " " + highWall.getX() + " " + SmoldynVCellMapper.SmoldynKeyword.p);
        } else {
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.low_wall + " 0 " + lowWall.getX() + " " + (compartSubDomain0.getBoundaryConditionXm().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.high_wall + " 0 " + highWall.getX() + " " + (compartSubDomain0.getBoundaryConditionXp().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
        }
        if (dimension > 1) {
            // y
            if (compartSubDomain0.getBoundaryConditionYm().isPERIODIC()) {
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 1 " + lowWall.getY() + " " + highWall.getY() + " " + SmoldynVCellMapper.SmoldynKeyword.p);
            } else {
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.low_wall + " 1 " + lowWall.getY() + " " + (compartSubDomain0.getBoundaryConditionYm().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.high_wall + " 1 " + highWall.getY() + " " + (compartSubDomain0.getBoundaryConditionYp().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
            }
            if (dimension > 2) {
                // z
                if (compartSubDomain0.getBoundaryConditionZm().isPERIODIC()) {
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 2 " + lowWall.getZ() + " " + highWall.getZ() + " " + SmoldynVCellMapper.SmoldynKeyword.p);
                } else {
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.low_wall + " 2 " + lowWall.getZ() + " " + (compartSubDomain0.getBoundaryConditionZm().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.high_wall + " 2 " + highWall.getZ() + " " + (compartSubDomain0.getBoundaryConditionZp().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
                }
            }
        }
        printWriter.println();
    } else {
        // x
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 0 " + lowWall.getX() + " " + highWall.getX());
        if (dimension > 1) {
            // y
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 1 " + lowWall.getY() + " " + highWall.getY());
            if (dimension > 2) {
                // z
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 2 " + lowWall.getZ() + " " + highWall.getZ());
            }
        }
        printWriter.println();
        // bounding walls as surfaces
        // have to find boundary condition type
        ISize sampleSize = simulation.getMeshSpecification().getSamplingSize();
        int numX = sampleSize.getX();
        int numY = dimension < 2 ? 1 : sampleSize.getY();
        int numZ = dimension < 3 ? 1 : sampleSize.getZ();
        if (dimension > 2) {
            int[] k_wall = new int[] { 0, numZ - 1 };
            for (int k = 0; k < k_wall.length; k++) {
                for (int j = 0; j < numY; j++) {
                    for (int i = 0; i < numX; i++) {
                        int volIndex = k_wall[k] * numX * numY + j * numX + i;
                        for (SubVolume sv : subVolumes) {
                            // gather all the points in all the regions
                            GeometricRegion[] geometricRegions = geometrySurfaceDescription.getGeometricRegions(sv);
                            RegionInfo[] regionInfos = geometrySurfaceDescription.getRegionImage().getRegionInfos();
                            for (GeometricRegion gr : geometricRegions) {
                                VolumeGeometricRegion vgr = (VolumeGeometricRegion) gr;
                                for (RegionInfo ri : regionInfos) {
                                    if (ri.getRegionIndex() == vgr.getRegionID() && ri.isIndexInRegion(volIndex)) {
                                        boundaryZSubVolumes.add(sv);
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        if (dimension > 1) {
            int[] j_wall = new int[] { 0, numY - 1 };
            for (int k = 0; k < numZ; k++) {
                for (int j = 0; j < j_wall.length; j++) {
                    for (int i = 0; i < numX; i++) {
                        int volIndex = k * numX * numY + j_wall[j] * numX + i;
                        for (SubVolume sv : subVolumes) {
                            // gather all the points in all the regions
                            GeometricRegion[] geometricRegions = geometrySurfaceDescription.getGeometricRegions(sv);
                            RegionInfo[] regionInfos = geometrySurfaceDescription.getRegionImage().getRegionInfos();
                            for (GeometricRegion gr : geometricRegions) {
                                VolumeGeometricRegion vgr = (VolumeGeometricRegion) gr;
                                for (RegionInfo ri : regionInfos) {
                                    if (ri.getRegionIndex() == vgr.getRegionID() && ri.isIndexInRegion(volIndex)) {
                                        boundaryYSubVolumes.add(sv);
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        int[] i_wall = new int[] { 0, numX - 1 };
        for (int k = 0; k < numZ; k++) {
            for (int j = 0; j < numY; j++) {
                for (int i = 0; i < i_wall.length; i++) {
                    int volIndex = k * numX * numY + j * numX + i_wall[i];
                    for (SubVolume sv : subVolumes) {
                        // gather all the points in all the regions
                        GeometricRegion[] geometricRegions = geometrySurfaceDescription.getGeometricRegions(sv);
                        RegionInfo[] regionInfos = geometrySurfaceDescription.getRegionImage().getRegionInfos();
                        for (GeometricRegion gr : geometricRegions) {
                            VolumeGeometricRegion vgr = (VolumeGeometricRegion) gr;
                            for (RegionInfo ri : regionInfos) {
                                if (ri.getRegionIndex() == vgr.getRegionID() && ri.isIndexInRegion(volIndex)) {
                                    boundaryXSubVolumes.add(sv);
                                }
                            }
                        }
                    }
                }
            }
        }
        Set<SubVolume> boundarySubVolumes = new HashSet<SubVolume>();
        boundarySubVolumes.addAll(boundaryXSubVolumes);
        boundarySubVolumes.addAll(boundaryYSubVolumes);
        boundarySubVolumes.addAll(boundaryZSubVolumes);
        BoundaryConditionType[] computedBct = new BoundaryConditionType[dimension * 2];
        String[] smoldynBct = new String[dimension * 2];
        String[] wallNames = new String[] { "Xm", "Xp", "Ym", "Yp", "Zm", "Zp" };
        if (boundarySubVolumes.size() >= 1) {
            for (SubVolume sv : boundarySubVolumes) {
                CompartmentSubDomain csd = (CompartmentSubDomain) mathDesc.getSubDomain(sv.getName());
                BoundaryConditionType[] bct = new BoundaryConditionType[] { csd.getBoundaryConditionXm(), csd.getBoundaryConditionXp(), csd.getBoundaryConditionYm(), csd.getBoundaryConditionYp(), csd.getBoundaryConditionZm(), csd.getBoundaryConditionZp() };
                if (computedBct[0] == null) {
                    System.arraycopy(bct, 0, computedBct, 0, dimension * 2);
                    for (int i = 0; i < dimension * 2; i++) {
                        if (computedBct[i].isPERIODIC()) {
                            throw new SolverException("Models with both surfaces and periodic boundary conditions are not supported yet.");
                        }
                        smoldynBct[i] = computedBct[i].isDIRICHLET() ? SmoldynVCellMapper.SmoldynKeyword.absorb.name() : SmoldynVCellMapper.SmoldynKeyword.reflect.name();
                    }
                } else {
                    for (int i = 0; i < dimension * 2; i++) {
                        if (!computedBct[i].compareEqual(bct[i])) {
                            throw new SolverException(wallNames[i] + " wall has different boundary conditions");
                        }
                    }
                }
            }
        }
        printWriter.println("# bounding wall surface");
        // X walls
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.start_surface + " " + VCellSmoldynKeyword.bounding_wall_surface_X);
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + "(" + SmoldynVCellMapper.SmoldynKeyword.up + ") " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.reflect);
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.front + " " + smoldynBct[0]);
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.back + " " + smoldynBct[1]);
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.color + " " + SmoldynVCellMapper.SmoldynKeyword.both + " 1 1 1");
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_panels + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " 2");
        // yz walls
        switch(dimension) {
            case 1:
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +0 " + lowWall.getX());
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -0 " + highWall.getX());
                break;
            case 2:
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +0 " + lowWall.getX() + " " + lowWall.getY() + " " + extent.getY());
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -0 " + highWall.getX() + " " + lowWall.getY() + " " + extent.getY());
                break;
            case 3:
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +0 " + lowWall.getX() + " " + lowWall.getY() + " " + lowWall.getZ() + " " + extent.getY() + " " + extent.getZ());
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -0 " + highWall.getX() + " " + lowWall.getY() + " " + lowWall.getZ() + " " + extent.getY() + " " + extent.getZ());
                break;
        }
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.end_surface);
        printWriter.println();
        if (dimension > 1) {
            // Y walls
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.start_surface + " " + VCellSmoldynKeyword.bounding_wall_surface_Y);
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + "(" + SmoldynVCellMapper.SmoldynKeyword.up + ") " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.reflect);
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.front + " " + smoldynBct[2]);
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.back + " " + smoldynBct[3]);
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.color + " " + SmoldynVCellMapper.SmoldynKeyword.both + " 1 1 1");
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_panels + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " 2");
            // xz walls
            switch(dimension) {
                case 2:
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +1 " + lowWall.getX() + " " + lowWall.getY() + " " + extent.getX());
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -1 " + lowWall.getX() + " " + highWall.getY() + " " + extent.getX());
                    break;
                case 3:
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +1 " + lowWall.getX() + " " + lowWall.getY() + " " + lowWall.getZ() + " " + extent.getX() + " " + extent.getZ());
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -1 " + lowWall.getX() + " " + highWall.getY() + " " + lowWall.getZ() + " " + extent.getX() + " " + extent.getZ());
                    break;
            }
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.end_surface);
            printWriter.println();
            if (dimension > 2) {
                // Z walls
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.start_surface + " " + VCellSmoldynKeyword.bounding_wall_surface_Z);
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + "(" + SmoldynVCellMapper.SmoldynKeyword.up + ") " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.reflect);
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.front + " " + smoldynBct[4]);
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.back + " " + smoldynBct[5]);
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.color + " " + SmoldynVCellMapper.SmoldynKeyword.both + " 1 1 1");
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_panels + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " 2");
                // xy walls
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +2 " + lowWall.getX() + " " + lowWall.getY() + " " + lowWall.getZ() + " " + extent.getX() + " " + extent.getY());
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -2 " + lowWall.getX() + " " + lowWall.getY() + " " + highWall.getZ() + " " + extent.getX() + " " + extent.getY());
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.end_surface);
                printWriter.println();
            }
        }
    }
}