Examples with BorderCellInfo org.vcell.vis.chombo.ChomboBoundaries.BorderCellInfo used on opensource projects

Example 1 with BorderCellInfo

use of org.vcell.vis.chombo.ChomboBoundaries.BorderCellInfo in project vcell by virtualcell.

the class ChomboMeshMapping method cropVoxels.

private void cropVoxels(VisMesh visMesh, ChomboBoundaries chomboBoundaries, ChomboCombinedVolumeMembraneDomain chomboCombinedVolumeMembraneDomain) {
    if (visMesh.getDimension() != 3) {
        throw new RuntimeException("expecting 3D mesh");
    }
    List<VisVoxel> origVoxelList = visMesh.getVisVoxels();
    ArrayList<VisVoxel> newVoxelList = new ArrayList<VisVoxel>();
    List<VisPoint> points = visMesh.getPoints();
    List<VisSurfaceTriangle> triangles = visMesh.getSurfaceTriangles();
    for (VisVoxel visVoxel : origVoxelList) {
        List<Integer> polyhedronPointIndices = visVoxel.getPointIndices();
        if (visVoxel.getChomboVolumeIndex().getFraction() < 1.0) {
            int p0 = polyhedronPointIndices.get(0);
            int p1 = polyhedronPointIndices.get(1);
            int p2 = polyhedronPointIndices.get(2);
            int p3 = polyhedronPointIndices.get(3);
            int p4 = polyhedronPointIndices.get(4);
            int p5 = polyhedronPointIndices.get(5);
            int p6 = polyhedronPointIndices.get(6);
            int p7 = polyhedronPointIndices.get(7);
            VisPoint vp0 = points.get(p0);
            VisPoint vp1 = points.get(p1);
            VisPoint vp2 = points.get(p2);
            VisPoint vp3 = points.get(p3);
            VisPoint vp4 = points.get(p4);
            VisPoint vp5 = points.get(p5);
            VisPoint vp6 = points.get(p6);
            VisPoint vp7 = points.get(p7);
            ArrayList<VisSurfaceTriangle> intersectingTriangles = new ArrayList<VisSurfaceTriangle>();
            for (VisSurfaceTriangle triangle : triangles) {
                boolean bInRange = true;
                for (int pi = 0; pi < 3; pi++) {
                    VisPoint tp = points.get(triangle.getPointIndices().get(pi));
                    if (!inLoHi(tp, vp0, vp7)) {
                        bInRange = false;
                    }
                }
                if (bInRange) {
                    intersectingTriangles.add(triangle);
                }
            }
            if (intersectingTriangles.size() == 0) {
                LG.info("fraction<1.0 but found no triangles");
                newVoxelList.add(visVoxel);
                continue;
            }
            // p6-------------------p7
            // /|                   /|
            // / |                  / |
            // p4-------------------p5  |
            // |  |                 |  |      face number         coordinates
            // |  |                 |  |
            // |  |                 |  |         5   3            z   y
            // |  p2................|..p3        |  /             |  /
            // | /                  | /          | /              | /
            // |/                   |/           |/               |/
            // p0-------------------p1       0 ---'---- 1          '----- x
            // /|
            // / |
            // 2  4
            BorderCellInfo borderCellInfo = chomboBoundaries.getMeshMetrics().getBorderCellInfo(intersectingTriangles.get(0).getChomboSurfaceIndex().getIndex());
            // 
            // have o flip the inside/outside if domain ordinal is > 0 ... note that "^" is the exclusive or ... to flip a bit
            // 
            VoxelPoint[] v = new VoxelPoint[8];
            for (int i = 0; i < 8; ++i) {
                int p = polyhedronPointIndices.get(i);
                VisPoint vp = points.get(p);
                v[i] = new VoxelPoint(p, vp, chomboCombinedVolumeMembraneDomain.shouldIncludeVertex(borderCellInfo.isVertexInPhase1(i)));
            }
            // choosing an arbitrary face (A,B,C,D) see below
            // 
            // pA   pB
            // 
            // pD   pC
            // 
            // face 0 (X-)
            VoxelFace face0 = new VoxelFace(Face.Xm, v[0], v[4], v[6], v[2]);
            // face 1 (X+)
            VoxelFace face1 = new VoxelFace(Face.Xp, v[1], v[3], v[7], v[5]);
            // face 2 (Y-)
            VoxelFace face2 = new VoxelFace(Face.Ym, v[0], v[1], v[5], v[4]);
            // face 3 (Y+)
            VoxelFace face3 = new VoxelFace(Face.Yp, v[2], v[3], v[7], v[6]);
            // face 4 (Z-)
            VoxelFace face4 = new VoxelFace(Face.Zm, v[0], v[2], v[3], v[1]);
            // face 5 (Z+)
            VoxelFace face5 = new VoxelFace(Face.Zp, v[4], v[5], v[7], v[6]);
            ClippedVoxel clippedVoxel = new ClippedVoxel(face0, face1, face2, face3, face4, face5);
            clippedVoxel.surfaceTriangles.addAll(intersectingTriangles);
            VisIrregularPolyhedron clippedPolyhedron = createClippedPolyhedron(clippedVoxel, visMesh, visVoxel);
            // VisIrregularPolyhedron clippedPolyhedron = new VisIrregularPolyhedron(visVoxel.getLevel(),visVoxel.getBoxNumber(),visVoxel.getBoxIndex(),visVoxel.getFraction());
            // clippedPolyhedron.addFace(new PolyhedronFace(new int[] { p0, p1, p4} ));
            // clippedPolyhedron.addFace(new PolyhedronFace(new int[] { p0, p2, p1} ));
            // clippedPolyhedron.addFace(new PolyhedronFace(new int[] { p0, p4, p2} ));
            // clippedPolyhedron.addFace(new PolyhedronFace(new int[] { p2, p4, p1} ));
            visMesh.addToIrregularPolyhedra(clippedPolyhedron);
        // VisTetrahedron[] delaunayTets = VtkGridUtils.createTetrahedra(clippedPolyhedron, visMesh);
        // for (VisTetrahedron tet : delaunayTets){
        // newPolyhedraList.add(tet);
        // }
        } else {
            // fraction >= 1.0
            newVoxelList.add(visVoxel);
        }
    }
    // for loop (orig polyhedra)
    visMesh.getVisVoxels().clear();
    visMesh.getVisVoxels().addAll(newVoxelList);
}

Example 2 with BorderCellInfo

use of org.vcell.vis.chombo.ChomboBoundaries.BorderCellInfo in project vcell by virtualcell.

the class ChomboFileReader method readMesh.

private static ChomboMeshData readMesh(String meshFileName, String vol0FileName) throws Exception {
    ChomboMesh chomboMesh = new ChomboMesh();
    if (H5.H5open() < 0) {
        throw new Exception("H5.H5open() failed");
    }
    FileFormat fileFormat = FileFormat.getFileFormat(FileFormat.FILE_TYPE_HDF5);
    if (fileFormat == null) {
        throw new Exception("FileFormat.getFileFormat(FileFormat.FILE_TYPE_HDF5) failed, returned null.");
    }
    FileFormat meshFile = fileFormat.createInstance(new File(meshFileName).getAbsolutePath(), FileFormat.READ);
    try {
        meshFile.open();
        DefaultMutableTreeNode meshRootNode = (DefaultMutableTreeNode) meshFile.getRootNode();
        Group meshRootGroup = (Group) meshRootNode.getUserObject();
        Group meshGroup = Hdf5Reader.getChildGroup(meshRootGroup, "mesh");
        chomboMesh.setDimension(Hdf5Reader.getIntAttribute(meshGroup, MESH_ATTR_DIMENSION));
        chomboMesh.setExtent(Hdf5Reader.getVect3DAttribute(meshGroup, MESH_ATTR_EXTENT, 1.0));
        chomboMesh.setOrigin(Hdf5Reader.getVect3DAttribute(meshGroup, MESH_ATTR_ORIGIN, 0.0));
        // it's very wasteful here, but what can I do?
        CartesianMeshChombo cartesianMeshChombo = CartesianMeshChombo.readMeshFile(new File(meshFileName));
        for (FeaturePhaseVol fpv : cartesianMeshChombo.getFeaturePhaseVols()) {
            chomboMesh.addFeaturePhase(fpv.feature, fpv.iphase);
        }
        // Hdf5Reader.DataColumn[] metricsColumns = Hdf5Reader.getDataTable(meshGroup,METRICS_DATASET);
        if (chomboMesh.getDimension() == 2) {
            Hdf5Reader.DataColumn[] segmentColumns = Hdf5Reader.getDataTable(meshGroup, "segments");
            Hdf5Reader.DataColumn[] verticesColumns = Hdf5Reader.getDataTable(meshGroup, "vertices");
            ChomboBoundaries boundaries = chomboMesh.getBoundaries();
            int numVertices = verticesColumns[0].getNumRows();
            int numSegments = segmentColumns[0].getNumRows();
            for (int i = 0; i < numVertices; i++) {
                double x = verticesColumns[0].getValue(i);
                double y = verticesColumns[1].getValue(i);
                double z = 0.0;
                boundaries.addPoint(new ChomboBoundaries.Point(x, y, z));
            }
            for (int i = 0; i < numSegments; i++) {
                int v1 = (int) segmentColumns[1].getValue(i);
                int v2 = (int) segmentColumns[2].getValue(i);
                // THIS COULD BE WRONG - is the chomboIndex one-to-one with the line segments? ... if not should be in the HDF5 file.
                int chomboIndex = i;
                boundaries.addSegment(new ChomboBoundaries.Segment(chomboIndex, v1, v2));
            }
        } else if (chomboMesh.getDimension() == 3) {
            Hdf5Reader.DataColumn[] surfaceTriangleColumns = Hdf5Reader.getDataTable(meshGroup, "surface triangles");
            ChomboBoundaries boundaries = chomboMesh.getBoundaries();
            int numTriangles = surfaceTriangleColumns[0].getNumRows();
            for (int row = 0; row < numTriangles; row++) {
                int index = (int) surfaceTriangleColumns[0].getValue(row);
                int faceNumber = (int) surfaceTriangleColumns[1].getValue(row);
                // not used currently
                int neighborIndex = (int) surfaceTriangleColumns[2].getValue(row);
                double x0 = surfaceTriangleColumns[3].getValue(row);
                double y0 = surfaceTriangleColumns[4].getValue(row);
                double z0 = surfaceTriangleColumns[5].getValue(row);
                int p0_index = boundaries.getOrCreatePoint(x0, y0, z0);
                double x1 = surfaceTriangleColumns[6].getValue(row);
                double y1 = surfaceTriangleColumns[7].getValue(row);
                double z1 = surfaceTriangleColumns[8].getValue(row);
                int p1_index = boundaries.getOrCreatePoint(x1, y1, z1);
                double x2 = surfaceTriangleColumns[9].getValue(row);
                double y2 = surfaceTriangleColumns[10].getValue(row);
                double z2 = surfaceTriangleColumns[11].getValue(row);
                int p2_index = boundaries.getOrCreatePoint(x2, y2, z2);
                Face face = Face.fromInteger(faceNumber);
                SurfaceTriangle surfaceTriangle = new SurfaceTriangle(index, face, p0_index, p1_index, p2_index);
                boundaries.addSurfaceTriangle(surfaceTriangle);
            }
            Hdf5Reader.DataColumn[] metricsColumns = Hdf5Reader.getDataTable(meshGroup, "membrane elements");
            MeshMetrics meshMetrics = boundaries.getMeshMetrics();
            int numMeshMetrics = metricsColumns[0].getNumRows();
            for (int row = 0; row < numMeshMetrics; row++) {
                int index = (int) metricsColumns[0].getValue(row);
                int level = (int) metricsColumns[1].getValue(row);
                // not used currently
                int i = (int) metricsColumns[2].getValue(row);
                // not used currently
                int j = (int) metricsColumns[3].getValue(row);
                // not used currently
                int k = (int) metricsColumns[4].getValue(row);
                double x = metricsColumns[5].getValue(row);
                double y = metricsColumns[6].getValue(row);
                double z = metricsColumns[7].getValue(row);
                Vect3D center = new Vect3D(x, y, z);
                double normalX = metricsColumns[8].getValue(row);
                double normalY = metricsColumns[9].getValue(row);
                double normalZ = metricsColumns[10].getValue(row);
                Vect3D normal = new Vect3D(normalX, normalY, normalZ);
                double volumeFraction = metricsColumns[11].getValue(row);
                double areaFraction = metricsColumns[12].getValue(row);
                int membraneId = (int) metricsColumns[13].getValue(row);
                int cornerPhaseMask = (int) metricsColumns[14].getValue(row);
                BorderCellInfo borderCellInfo = new BorderCellInfo(index, level, i, j, k, center, normal, volumeFraction, areaFraction, membraneId, cornerPhaseMask);
                meshMetrics.addBorderCellInfo(borderCellInfo);
            }
        } else {
            throw new Exception("failed to read chombo file, unexpected mesh dimension " + chomboMesh.getDimension());
        }
    } finally {
        meshFile.close();
    }
    FileFormat vol0File = fileFormat.createInstance(new File(vol0FileName).getAbsolutePath(), FileFormat.READ);
    try {
        vol0File.open();
        DefaultMutableTreeNode rootNode = (DefaultMutableTreeNode) vol0File.getRootNode();
        Group rootGroup = (Group) rootNode.getUserObject();
        Group level0Group = Hdf5Reader.getChildGroup(rootGroup, "level_0");
        double time = Hdf5Reader.getDoubleAttribute(level0Group, "time");
        ChomboMeshData chomboMeshData = new ChomboMeshData(chomboMesh, time);
        int numComponents = Hdf5Reader.getIntAttribute(rootGroup, "num_components");
        int numLevels = Hdf5Reader.getIntAttribute(rootGroup, "num_levels");
        int fractionComponentIndex = -1;
        for (int i = 0; i < numComponents; i++) {
            String compName = Hdf5Reader.getStringAttribute(rootGroup, "component_" + i);
            chomboMeshData.addComponentName(compName);
            if (compName.equals("fraction-0")) {
                fractionComponentIndex = i;
                chomboMeshData.setFraction0ComponentIndex(fractionComponentIndex);
            }
        }
        for (int i = 0; i < numLevels; i++) {
            Group levelGroup = Hdf5Reader.getChildGroup(rootGroup, "level_" + i);
            int refinement = 2;
            if (i == 0) {
                refinement = 1;
            }
            ChomboLevel chomboLevel = new ChomboLevel(chomboMesh, i, refinement);
            Hdf5Reader.DataColumn[] boxColumns = Hdf5Reader.getDataTable(levelGroup, "boxes");
            int[] lo_i, lo_j, lo_k, hi_i, hi_j, hi_k;
            if (chomboMesh.getDimension() == 2) {
                lo_i = ((Hdf5Reader.IntColumn) boxColumns[0]).data;
                lo_j = ((Hdf5Reader.IntColumn) boxColumns[1]).data;
                hi_i = ((Hdf5Reader.IntColumn) boxColumns[2]).data;
                hi_j = ((Hdf5Reader.IntColumn) boxColumns[3]).data;
                lo_k = new int[boxColumns[0].getNumRows()];
                hi_k = new int[boxColumns[0].getNumRows()];
            } else {
                lo_i = ((Hdf5Reader.IntColumn) boxColumns[0]).data;
                lo_j = ((Hdf5Reader.IntColumn) boxColumns[1]).data;
                lo_k = ((Hdf5Reader.IntColumn) boxColumns[2]).data;
                hi_i = ((Hdf5Reader.IntColumn) boxColumns[3]).data;
                hi_j = ((Hdf5Reader.IntColumn) boxColumns[4]).data;
                hi_k = ((Hdf5Reader.IntColumn) boxColumns[5]).data;
            }
            for (int b = 0; b < boxColumns[0].getNumRows(); b++) {
                ChomboBox chomboBox = new ChomboBox(chomboLevel, lo_i[b], hi_i[b], lo_j[b], hi_j[b], lo_k[b], hi_k[b], chomboMesh.getDimension());
                chomboLevel.addBox(chomboBox);
            }
            chomboMesh.addLevel(chomboLevel);
            // 
            // read the variables
            // 
            Hdf5Reader.DataColumn[] data = Hdf5Reader.getDataTable(levelGroup, "data:datatype=0");
            Hdf5Reader.DataColumn[] offsets = Hdf5Reader.getDataTable(levelGroup, "data:offsets=0");
            ChomboLevelData chomboLevelData = new ChomboLevelData(i, fractionComponentIndex, ((Hdf5Reader.DoubleColumn) data[0]).data, ((Hdf5Reader.LongColumn) offsets[0]).data);
            chomboMeshData.addLevelData(chomboLevelData);
        }
        readMembraneVarData(chomboMeshData, rootGroup);
        return chomboMeshData;
    } finally {
        vol0File.close();
    }
}