Search in sources :

Example 1 with AnalyticSubVolume

use of cbit.vcell.geometry.AnalyticSubVolume in project vcell by virtualcell.

the class SBMLImporter method addGeometry.

protected void addGeometry() {
    // get a Geometry object via SpatialModelPlugin object.
    org.sbml.jsbml.ext.spatial.Geometry sbmlGeometry = getSbmlGeometry();
    if (sbmlGeometry == null) {
        return;
    }
    int dimension = 0;
    Origin vcOrigin = null;
    Extent vcExtent = null;
    {
        // local code block
        // get a CoordComponent object via the Geometry object.
        ListOf<CoordinateComponent> listOfCoordComps = sbmlGeometry.getListOfCoordinateComponents();
        if (listOfCoordComps == null) {
            throw new RuntimeException("Cannot have 0 coordinate compartments in geometry");
        }
        // coord component
        double ox = 0.0;
        double oy = 0.0;
        double oz = 0.0;
        double ex = 1.0;
        double ey = 1.0;
        double ez = 1.0;
        for (CoordinateComponent coordComponent : listOfCoordComps) {
            double minValue = coordComponent.getBoundaryMinimum().getValue();
            double maxValue = coordComponent.getBoundaryMaximum().getValue();
            switch(coordComponent.getType()) {
                case cartesianX:
                    {
                        ox = minValue;
                        ex = maxValue - minValue;
                        break;
                    }
                case cartesianY:
                    {
                        oy = minValue;
                        ey = maxValue - minValue;
                        break;
                    }
                case cartesianZ:
                    {
                        oz = minValue;
                        ez = maxValue - minValue;
                        break;
                    }
            }
            dimension++;
        }
        vcOrigin = new Origin(ox, oy, oz);
        vcExtent = new Extent(ex, ey, ez);
    }
    // from geometry definition, find out which type of geometry : image or
    // analytic or CSG
    AnalyticGeometry analyticGeometryDefinition = null;
    CSGeometry csGeometry = null;
    SampledFieldGeometry segmentedSampledFieldGeometry = null;
    SampledFieldGeometry distanceMapSampledFieldGeometry = null;
    ParametricGeometry parametricGeometry = null;
    for (int i = 0; i < sbmlGeometry.getListOfGeometryDefinitions().size(); i++) {
        GeometryDefinition gd_temp = sbmlGeometry.getListOfGeometryDefinitions().get(i);
        if (!gd_temp.isSetIsActive()) {
            continue;
        }
        if (gd_temp instanceof AnalyticGeometry) {
            analyticGeometryDefinition = (AnalyticGeometry) gd_temp;
        } else if (gd_temp instanceof SampledFieldGeometry) {
            SampledFieldGeometry sfg = (SampledFieldGeometry) gd_temp;
            String sfn = sfg.getSampledField();
            ListOf<SampledField> sampledFields = sbmlGeometry.getListOfSampledFields();
            if (sampledFields.size() > 1) {
                throw new RuntimeException("only one sampled field supported");
            }
            InterpolationKind ik = sampledFields.get(0).getInterpolationType();
            switch(ik) {
                case linear:
                    distanceMapSampledFieldGeometry = sfg;
                    break;
                case nearestneighbor:
                    segmentedSampledFieldGeometry = sfg;
                    break;
                default:
                    lg.warn("Unsupported " + sampledFields.get(0).getName() + " interpolation type " + ik);
            }
        } else if (gd_temp instanceof CSGeometry) {
            csGeometry = (CSGeometry) gd_temp;
        } else if (gd_temp instanceof ParametricGeometry) {
            parametricGeometry = (ParametricGeometry) gd_temp;
        } else {
            throw new RuntimeException("unsupported geometry definition type " + gd_temp.getClass().getSimpleName());
        }
    }
    if (analyticGeometryDefinition == null && segmentedSampledFieldGeometry == null && distanceMapSampledFieldGeometry == null && csGeometry == null) {
        throw new SBMLImportException("VCell supports only Analytic, Image based (segmentd or distance map) or Constructed Solid Geometry at this time.");
    }
    GeometryDefinition selectedGeometryDefinition = null;
    if (csGeometry != null) {
        selectedGeometryDefinition = csGeometry;
    } else if (analyticGeometryDefinition != null) {
        selectedGeometryDefinition = analyticGeometryDefinition;
    } else if (segmentedSampledFieldGeometry != null) {
        selectedGeometryDefinition = segmentedSampledFieldGeometry;
    } else if (distanceMapSampledFieldGeometry != null) {
        selectedGeometryDefinition = distanceMapSampledFieldGeometry;
    } else if (parametricGeometry != null) {
        selectedGeometryDefinition = parametricGeometry;
    } else {
        throw new SBMLImportException("no geometry definition found");
    }
    Geometry vcGeometry = null;
    if (selectedGeometryDefinition == analyticGeometryDefinition || selectedGeometryDefinition == csGeometry) {
        vcGeometry = new Geometry("spatialGeom", dimension);
    } else if (selectedGeometryDefinition == distanceMapSampledFieldGeometry || selectedGeometryDefinition == segmentedSampledFieldGeometry) {
        SampledFieldGeometry sfg = (SampledFieldGeometry) selectedGeometryDefinition;
        // get image from sampledFieldGeometry
        // get a sampledVol object via the listOfSampledVol (from
        // SampledGeometry) object.
        // gcw gcw gcw
        String sfn = sfg.getSampledField();
        SampledField sf = null;
        for (SampledField sampledField : sbmlGeometry.getListOfSampledFields()) {
            if (sampledField.getSpatialId().equals(sfn)) {
                sf = sampledField;
            }
        }
        int numX = sf.getNumSamples1();
        int numY = sf.getNumSamples2();
        int numZ = sf.getNumSamples3();
        int[] samples = new int[sf.getSamplesLength()];
        StringTokenizer tokens = new StringTokenizer(sf.getSamples(), " ");
        int count = 0;
        while (tokens.hasMoreTokens()) {
            int sample = Integer.parseInt(tokens.nextToken());
            samples[count++] = sample;
        }
        byte[] imageInBytes = new byte[samples.length];
        if (selectedGeometryDefinition == distanceMapSampledFieldGeometry) {
            // 
            for (int i = 0; i < imageInBytes.length; i++) {
                // if (interpolation(samples[i])<0){
                if (samples[i] < 0) {
                    imageInBytes[i] = -1;
                } else {
                    imageInBytes[i] = 1;
                }
            }
        } else {
            for (int i = 0; i < imageInBytes.length; i++) {
                imageInBytes[i] = (byte) samples[i];
            }
        }
        try {
            // System.out.println("ident " + sf.getId() + " " + sf.getName());
            VCImage vcImage = null;
            CompressionKind ck = sf.getCompression();
            DataKind dk = sf.getDataType();
            if (ck == CompressionKind.deflated) {
                vcImage = new VCImageCompressed(null, imageInBytes, vcExtent, numX, numY, numZ);
            } else {
                switch(dk) {
                    case UINT8:
                    case UINT16:
                    case UINT32:
                        vcImage = new VCImageUncompressed(null, imageInBytes, vcExtent, numX, numY, numZ);
                    default:
                }
            }
            if (vcImage == null) {
                throw new SbmlException("Unsupported type combination " + ck + ", " + dk + " for sampled field " + sf.getName());
            }
            vcImage.setName(sf.getId());
            ListOf<SampledVolume> sampledVolumes = sfg.getListOfSampledVolumes();
            final int numSampledVols = sampledVolumes.size();
            if (numSampledVols == 0) {
                throw new RuntimeException("Cannot have 0 sampled volumes in sampledField (image_based) geometry");
            }
            // check to see if values are uniquely integer , add set up scaling if necessary
            double scaleFactor = checkPixelScaling(sampledVolumes, 1);
            if (scaleFactor != 1) {
                double checkScaleFactor = checkPixelScaling(sampledVolumes, scaleFactor);
                VCAssert.assertTrue(checkScaleFactor != scaleFactor, "Scale factor check failed");
            }
            VCPixelClass[] vcpixelClasses = new VCPixelClass[numSampledVols];
            // get pixel classes for geometry
            for (int i = 0; i < numSampledVols; i++) {
                SampledVolume sVol = sampledVolumes.get(i);
                // from subVolume, get pixelClass?
                final int scaled = (int) (scaleFactor * sVol.getSampledValue());
                vcpixelClasses[i] = new VCPixelClass(null, sVol.getDomainType(), scaled);
            }
            vcImage.setPixelClasses(vcpixelClasses);
            // now create image geometry
            vcGeometry = new Geometry("spatialGeom", vcImage);
        } catch (Exception e) {
            e.printStackTrace(System.out);
            throw new RuntimeException("Unable to create image from SampledFieldGeometry : " + e.getMessage());
        }
    }
    GeometrySpec vcGeometrySpec = vcGeometry.getGeometrySpec();
    vcGeometrySpec.setOrigin(vcOrigin);
    try {
        vcGeometrySpec.setExtent(vcExtent);
    } catch (PropertyVetoException e) {
        e.printStackTrace(System.out);
        throw new SBMLImportException("Unable to set extent on VC geometry : " + e.getMessage(), e);
    }
    // get listOfDomainTypes via the Geometry object.
    ListOf<DomainType> listOfDomainTypes = sbmlGeometry.getListOfDomainTypes();
    if (listOfDomainTypes == null || listOfDomainTypes.size() < 1) {
        throw new SBMLImportException("Cannot have 0 domainTypes in geometry");
    }
    // get a listOfDomains via the Geometry object.
    ListOf<Domain> listOfDomains = sbmlGeometry.getListOfDomains();
    if (listOfDomains == null || listOfDomains.size() < 1) {
        throw new SBMLImportException("Cannot have 0 domains in geometry");
    }
    // ListOfGeometryDefinitions listOfGeomDefns =
    // sbmlGeometry.getListOfGeometryDefinitions();
    // if ((listOfGeomDefns == null) ||
    // (sbmlGeometry.getNumGeometryDefinitions() > 1)) {
    // throw new
    // RuntimeException("Can have only 1 geometry definition in geometry");
    // }
    // use the boolean bAnalytic to create the right kind of subvolume.
    // First match the somVol=domainTypes for spDim=3. Deal witl spDim=2
    // afterwards.
    GeometrySurfaceDescription vcGsd = vcGeometry.getGeometrySurfaceDescription();
    Vector<DomainType> surfaceClassDomainTypesVector = new Vector<DomainType>();
    try {
        for (DomainType dt : listOfDomainTypes) {
            if (dt.getSpatialDimensions() == 3) {
                // subvolume
                if (selectedGeometryDefinition == analyticGeometryDefinition) {
                    // will set expression later - when reading in Analytic
                    // Volumes in GeometryDefinition
                    vcGeometrySpec.addSubVolume(new AnalyticSubVolume(dt.getId(), new Expression(1.0)));
                } else {
                // add SubVolumes later for CSG and Image-based
                }
            } else if (dt.getSpatialDimensions() == 2) {
                surfaceClassDomainTypesVector.add(dt);
            }
        }
        // analytic vol is needed to get the expression for subVols
        if (selectedGeometryDefinition == analyticGeometryDefinition) {
            // get an analyticVol object via the listOfAnalyticVol (from
            // AnalyticGeometry) object.
            ListOf<AnalyticVolume> aVolumes = analyticGeometryDefinition.getListOfAnalyticVolumes();
            if (aVolumes.size() < 1) {
                throw new SBMLImportException("Cannot have 0 Analytic volumes in analytic geometry");
            }
            for (AnalyticVolume analyticVol : aVolumes) {
                // get subVol from VC geometry using analyticVol spatialId;
                // set its expr using analyticVol's math.
                SubVolume vcSubvolume = vcGeometrySpec.getSubVolume(analyticVol.getDomainType());
                CastInfo<AnalyticSubVolume> ci = BeanUtils.attemptCast(AnalyticSubVolume.class, vcSubvolume);
                if (!ci.isGood()) {
                    throw new RuntimeException("analytic volume '" + analyticVol.getId() + "' does not map to any VC subvolume.");
                }
                AnalyticSubVolume asv = ci.get();
                try {
                    Expression subVolExpr = getExpressionFromFormula(analyticVol.getMath());
                    asv.setExpression(subVolExpr);
                } catch (ExpressionException e) {
                    e.printStackTrace(System.out);
                    throw new SBMLImportException("Unable to set expression on subVolume '" + asv.getName() + "'. " + e.getMessage(), e);
                }
            }
        }
        SampledFieldGeometry sfg = BeanUtils.downcast(SampledFieldGeometry.class, selectedGeometryDefinition);
        if (sfg != null) {
            ListOf<SampledVolume> sampledVolumes = sfg.getListOfSampledVolumes();
            int numSampledVols = sampledVolumes.size();
            if (numSampledVols == 0) {
                throw new SBMLImportException("Cannot have 0 sampled volumes in sampledField (image_based) geometry");
            }
            VCPixelClass[] vcpixelClasses = new VCPixelClass[numSampledVols];
            ImageSubVolume[] vcImageSubVols = new ImageSubVolume[numSampledVols];
            // get pixel classes for geometry
            int idx = 0;
            for (SampledVolume sVol : sampledVolumes) {
                // from subVolume, get pixelClass?
                final String name = sVol.getDomainType();
                final int pixelValue = SBMLUtils.ignoreZeroFraction(sVol.getSampledValue());
                VCPixelClass pc = new VCPixelClass(null, name, pixelValue);
                vcpixelClasses[idx] = pc;
                // Create the new Image SubVolume - use index of this for
                // loop as 'handle' for ImageSubVol?
                ImageSubVolume isv = new ImageSubVolume(null, pc, idx);
                isv.setName(name);
                vcImageSubVols[idx++] = isv;
            }
            vcGeometry.getGeometrySpec().setSubVolumes(vcImageSubVols);
        }
        if (selectedGeometryDefinition == csGeometry) {
            ListOf<org.sbml.jsbml.ext.spatial.CSGObject> listOfcsgObjs = csGeometry.getListOfCSGObjects();
            ArrayList<org.sbml.jsbml.ext.spatial.CSGObject> sbmlCSGs = new ArrayList<org.sbml.jsbml.ext.spatial.CSGObject>(listOfcsgObjs);
            // we want the CSGObj with highest ordinal to be the first
            // element in the CSG subvols array.
            Collections.sort(sbmlCSGs, new Comparator<org.sbml.jsbml.ext.spatial.CSGObject>() {

                @Override
                public int compare(org.sbml.jsbml.ext.spatial.CSGObject lhs, org.sbml.jsbml.ext.spatial.CSGObject rhs) {
                    // minus one to reverse sort
                    return -1 * Integer.compare(lhs.getOrdinal(), rhs.getOrdinal());
                }
            });
            int n = sbmlCSGs.size();
            CSGObject[] vcCSGSubVolumes = new CSGObject[n];
            for (int i = 0; i < n; i++) {
                org.sbml.jsbml.ext.spatial.CSGObject sbmlCSGObject = sbmlCSGs.get(i);
                CSGObject vcellCSGObject = new CSGObject(null, sbmlCSGObject.getDomainType(), i);
                vcellCSGObject.setRoot(getVCellCSGNode(sbmlCSGObject.getCSGNode()));
            }
            vcGeometry.getGeometrySpec().setSubVolumes(vcCSGSubVolumes);
        }
        // Call geom.geomSurfDesc.updateAll() to automatically generate
        // surface classes.
        // vcGsd.updateAll();
        vcGeometry.precomputeAll(new GeometryThumbnailImageFactoryAWT(), true, true);
    } catch (Exception e) {
        e.printStackTrace(System.out);
        throw new SBMLImportException("Unable to create VC subVolumes from SBML domainTypes : " + e.getMessage(), e);
    }
    // should now map each SBML domain to right VC geometric region.
    GeometricRegion[] vcGeomRegions = vcGsd.getGeometricRegions();
    ISize sampleSize = vcGsd.getVolumeSampleSize();
    RegionInfo[] regionInfos = vcGsd.getRegionImage().getRegionInfos();
    int numX = sampleSize.getX();
    int numY = sampleSize.getY();
    int numZ = sampleSize.getZ();
    double ox = vcOrigin.getX();
    double oy = vcOrigin.getY();
    double oz = vcOrigin.getZ();
    for (Domain domain : listOfDomains) {
        String domainType = domain.getDomainType();
        InteriorPoint interiorPt = domain.getListOfInteriorPoints().get(0);
        if (interiorPt == null) {
            DomainType currDomainType = null;
            for (DomainType dt : sbmlGeometry.getListOfDomainTypes()) {
                if (dt.getSpatialId().equals(domainType)) {
                    currDomainType = dt;
                }
            }
            if (currDomainType.getSpatialDimensions() == 2) {
                continue;
            }
        }
        Coordinate sbmlInteriorPtCoord = new Coordinate(interiorPt.getCoord1(), interiorPt.getCoord2(), interiorPt.getCoord3());
        for (int j = 0; j < vcGeomRegions.length; j++) {
            if (vcGeomRegions[j] instanceof VolumeGeometricRegion) {
                int regionID = ((VolumeGeometricRegion) vcGeomRegions[j]).getRegionID();
                for (int k = 0; k < regionInfos.length; k++) {
                    // (using gemoRegion regionID).
                    if (regionInfos[k].getRegionIndex() == regionID) {
                        int volIndx = 0;
                        Coordinate nearestPtCoord = null;
                        double minDistance = Double.MAX_VALUE;
                        // represented by SBML 'domain[i]'.
                        for (int z = 0; z < numZ; z++) {
                            for (int y = 0; y < numY; y++) {
                                for (int x = 0; x < numX; x++) {
                                    if (regionInfos[k].isIndexInRegion(volIndx)) {
                                        double unit_z = (numZ > 1) ? ((double) z) / (numZ - 1) : 0.5;
                                        double coordZ = oz + vcExtent.getZ() * unit_z;
                                        double unit_y = (numY > 1) ? ((double) y) / (numY - 1) : 0.5;
                                        double coordY = oy + vcExtent.getY() * unit_y;
                                        double unit_x = (numX > 1) ? ((double) x) / (numX - 1) : 0.5;
                                        double coordX = ox + vcExtent.getX() * unit_x;
                                        // for now, find the shortest dist
                                        // coord. Can refine algo later.
                                        Coordinate vcCoord = new Coordinate(coordX, coordY, coordZ);
                                        double distance = sbmlInteriorPtCoord.distanceTo(vcCoord);
                                        if (distance < minDistance) {
                                            minDistance = distance;
                                            nearestPtCoord = vcCoord;
                                        }
                                    }
                                    volIndx++;
                                }
                            // end - for x
                            }
                        // end - for y
                        }
                        // with domain name
                        if (nearestPtCoord != null) {
                            GeometryClass geomClassSBML = vcGeometry.getGeometryClass(domainType);
                            // we know vcGeometryReg[j] is a VolGeomRegion
                            GeometryClass geomClassVC = ((VolumeGeometricRegion) vcGeomRegions[j]).getSubVolume();
                            if (geomClassSBML.compareEqual(geomClassVC)) {
                                vcGeomRegions[j].setName(domain.getId());
                            }
                        }
                    }
                // end if (regInfoIndx = regId)
                }
            // end - for regInfo
            }
        }
    // end for - vcGeomRegions
    }
    // deal with surfaceClass:spDim2-domainTypes
    for (int i = 0; i < surfaceClassDomainTypesVector.size(); i++) {
        DomainType surfaceClassDomainType = surfaceClassDomainTypesVector.elementAt(i);
        // 'surfaceClassDomainType'
        for (Domain d : listOfDomains) {
            if (d.getDomainType().equals(surfaceClassDomainType.getId())) {
                // get the adjacent domains of this 'surface' domain
                // (surface domain + its 2 adj vol domains)
                Set<Domain> adjacentDomainsSet = getAssociatedAdjacentDomains(sbmlGeometry, d);
                // get the domain types of the adjacent domains in SBML and
                // store the corresponding subVol counterparts from VC for
                // adj vol domains
                Vector<SubVolume> adjacentSubVolumesVector = new Vector<SubVolume>();
                Vector<VolumeGeometricRegion> adjVolGeomRegionsVector = new Vector<VolumeGeometricRegion>();
                Iterator<Domain> iterator = adjacentDomainsSet.iterator();
                while (iterator.hasNext()) {
                    Domain dom = iterator.next();
                    DomainType dt = getBySpatialID(sbmlGeometry.getListOfDomainTypes(), dom.getDomainType());
                    if (dt.getSpatialDimensions() == 3) {
                        // for domain type with sp. dim = 3, get
                        // correspoinding subVol from VC geometry.
                        GeometryClass gc = vcGeometry.getGeometryClass(dt.getId());
                        adjacentSubVolumesVector.add((SubVolume) gc);
                        // store volGeomRegions corresponding to this (vol)
                        // geomClass in adjVolGeomRegionsVector : this
                        // should return ONLY 1 region for subVol.
                        GeometricRegion[] geomRegion = vcGsd.getGeometricRegions(gc);
                        adjVolGeomRegionsVector.add((VolumeGeometricRegion) geomRegion[0]);
                    }
                }
                // there should be only 2 subVols in this vector
                if (adjacentSubVolumesVector.size() != 2) {
                    throw new RuntimeException("Cannot have more or less than 2 subvolumes that are adjacent to surface (membrane) '" + d.getId() + "'");
                }
                // get the surface class with these 2 adj subVols. Set its
                // name to that of 'surfaceClassDomainType'
                SurfaceClass surfacClass = vcGsd.getSurfaceClass(adjacentSubVolumesVector.get(0), adjacentSubVolumesVector.get(1));
                surfacClass.setName(surfaceClassDomainType.getSpatialId());
                // get surfaceGeometricRegion that has adjVolGeomRegions as
                // its adjacent vol geom regions and set its name from
                // domain 'd'
                SurfaceGeometricRegion surfaceGeomRegion = getAssociatedSurfaceGeometricRegion(vcGsd, adjVolGeomRegionsVector);
                if (surfaceGeomRegion != null) {
                    surfaceGeomRegion.setName(d.getId());
                }
            }
        // end if - domain.domainType == surfaceClassDomainType
        }
    // end for - numDomains
    }
    // structureMappings in VC from compartmentMappings in SBML
    try {
        // set geometry first and then set structureMappings?
        vcBioModel.getSimulationContext(0).setGeometry(vcGeometry);
        // update simContextName ...
        vcBioModel.getSimulationContext(0).setName(vcBioModel.getSimulationContext(0).getName() + "_" + vcGeometry.getName());
        Model vcModel = vcBioModel.getSimulationContext(0).getModel();
        ModelUnitSystem vcModelUnitSystem = vcModel.getUnitSystem();
        Vector<StructureMapping> structMappingsVector = new Vector<StructureMapping>();
        SpatialCompartmentPlugin cplugin = null;
        for (int i = 0; i < sbmlModel.getNumCompartments(); i++) {
            Compartment c = sbmlModel.getCompartment(i);
            String cname = c.getName();
            cplugin = (SpatialCompartmentPlugin) c.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
            CompartmentMapping compMapping = cplugin.getCompartmentMapping();
            if (compMapping != null) {
                // final String id = compMapping.getId();
                // final String name = compMapping.getName();
                CastInfo<Structure> ci = SBMLHelper.getTypedStructure(Structure.class, vcModel, cname);
                if (ci.isGood()) {
                    Structure struct = ci.get();
                    String domainType = compMapping.getDomainType();
                    GeometryClass geometryClass = vcGeometry.getGeometryClass(domainType);
                    double unitSize = compMapping.getUnitSize();
                    Feature feat = BeanUtils.downcast(Feature.class, struct);
                    if (feat != null) {
                        FeatureMapping featureMapping = new FeatureMapping(feat, vcBioModel.getSimulationContext(0), vcModelUnitSystem);
                        featureMapping.setGeometryClass(geometryClass);
                        if (geometryClass instanceof SubVolume) {
                            featureMapping.getVolumePerUnitVolumeParameter().setExpression(new Expression(unitSize));
                        } else if (geometryClass instanceof SurfaceClass) {
                            featureMapping.getVolumePerUnitAreaParameter().setExpression(new Expression(unitSize));
                        }
                        structMappingsVector.add(featureMapping);
                    } else if (struct instanceof Membrane) {
                        MembraneMapping membraneMapping = new MembraneMapping((Membrane) struct, vcBioModel.getSimulationContext(0), vcModelUnitSystem);
                        membraneMapping.setGeometryClass(geometryClass);
                        if (geometryClass instanceof SubVolume) {
                            membraneMapping.getAreaPerUnitVolumeParameter().setExpression(new Expression(unitSize));
                        } else if (geometryClass instanceof SurfaceClass) {
                            membraneMapping.getAreaPerUnitAreaParameter().setExpression(new Expression(unitSize));
                        }
                        structMappingsVector.add(membraneMapping);
                    }
                }
            }
        }
        StructureMapping[] structMappings = structMappingsVector.toArray(new StructureMapping[0]);
        vcBioModel.getSimulationContext(0).getGeometryContext().setStructureMappings(structMappings);
        // if type from SBML parameter Boundary Condn is not the same as the
        // boundary type of the
        // structureMapping of structure of paramSpContext, set the boundary
        // condn type of the structureMapping
        // to the value of 'type' from SBML parameter Boundary Condn.
        ListOf<Parameter> listOfGlobalParams = sbmlModel.getListOfParameters();
        for (Parameter sbmlGlobalParam : sbmlModel.getListOfParameters()) {
            SpatialParameterPlugin spplugin = (SpatialParameterPlugin) sbmlGlobalParam.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
            ParameterType paramType = spplugin.getParamType();
            if (!(paramType instanceof BoundaryCondition)) {
                continue;
            }
            BoundaryCondition bCondn = (BoundaryCondition) paramType;
            if (bCondn.isSetVariable()) {
                // get the var of boundaryCondn; find appropriate spContext
                // in vcell;
                SpeciesContext paramSpContext = vcBioModel.getSimulationContext(0).getModel().getSpeciesContext(bCondn.getVariable());
                if (paramSpContext != null) {
                    Structure s = paramSpContext.getStructure();
                    StructureMapping sm = vcBioModel.getSimulationContext(0).getGeometryContext().getStructureMapping(s);
                    if (sm != null) {
                        BoundaryConditionType bct = null;
                        switch(bCondn.getType()) {
                            case Dirichlet:
                                {
                                    bct = BoundaryConditionType.DIRICHLET;
                                    break;
                                }
                            case Neumann:
                                {
                                    bct = BoundaryConditionType.NEUMANN;
                                    break;
                                }
                            case Robin_inwardNormalGradientCoefficient:
                            case Robin_sum:
                            case Robin_valueCoefficient:
                            default:
                                throw new RuntimeException("boundary condition type " + bCondn.getType().name() + " not supported");
                        }
                        for (CoordinateComponent coordComp : getSbmlGeometry().getListOfCoordinateComponents()) {
                            if (bCondn.getSpatialRef().equals(coordComp.getBoundaryMinimum().getSpatialId())) {
                                switch(coordComp.getType()) {
                                    case cartesianX:
                                        {
                                            sm.setBoundaryConditionTypeXm(bct);
                                        }
                                    case cartesianY:
                                        {
                                            sm.setBoundaryConditionTypeYm(bct);
                                        }
                                    case cartesianZ:
                                        {
                                            sm.setBoundaryConditionTypeZm(bct);
                                        }
                                }
                            }
                            if (bCondn.getSpatialRef().equals(coordComp.getBoundaryMaximum().getSpatialId())) {
                                switch(coordComp.getType()) {
                                    case cartesianX:
                                        {
                                            sm.setBoundaryConditionTypeXm(bct);
                                        }
                                    case cartesianY:
                                        {
                                            sm.setBoundaryConditionTypeYm(bct);
                                        }
                                    case cartesianZ:
                                        {
                                            sm.setBoundaryConditionTypeZm(bct);
                                        }
                                }
                            }
                        }
                    } else // sm != null
                    {
                        logger.sendMessage(VCLogger.Priority.MediumPriority, VCLogger.ErrorType.OverallWarning, "No structure " + s.getName() + " requested by species context " + paramSpContext.getName());
                    }
                }
            // end if (paramSpContext != null)
            }
        // end if (bCondn.isSetVar())
        }
        // end for (sbmlModel.numParams)
        vcBioModel.getSimulationContext(0).getGeometryContext().refreshStructureMappings();
        vcBioModel.getSimulationContext(0).refreshSpatialObjects();
    } catch (Exception e) {
        e.printStackTrace(System.out);
        throw new SBMLImportException("Unable to create VC structureMappings from SBML compartment mappings : " + e.getMessage(), e);
    }
}
Also used : Origin(org.vcell.util.Origin) VCPixelClass(cbit.image.VCPixelClass) MembraneMapping(cbit.vcell.mapping.MembraneMapping) DataKind(org.sbml.jsbml.ext.spatial.DataKind) ArrayList(java.util.ArrayList) BoundaryConditionType(cbit.vcell.math.BoundaryConditionType) SpeciesContext(cbit.vcell.model.SpeciesContext) Feature(cbit.vcell.model.Feature) GeometryDefinition(org.sbml.jsbml.ext.spatial.GeometryDefinition) SubVolume(cbit.vcell.geometry.SubVolume) ImageSubVolume(cbit.vcell.geometry.ImageSubVolume) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) Vector(java.util.Vector) CoordinateComponent(org.sbml.jsbml.ext.spatial.CoordinateComponent) SimulationContext(cbit.vcell.mapping.SimulationContext) SpeciesContext(cbit.vcell.model.SpeciesContext) IssueContext(org.vcell.util.IssueContext) ReactionContext(cbit.vcell.mapping.ReactionContext) CompressionKind(org.sbml.jsbml.ext.spatial.CompressionKind) VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion) PropertyVetoException(java.beans.PropertyVetoException) ModelPropertyVetoException(cbit.vcell.model.ModelPropertyVetoException) Coordinate(org.vcell.util.Coordinate) BoundaryCondition(org.sbml.jsbml.ext.spatial.BoundaryCondition) SbmlException(org.vcell.sbml.SbmlException) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) SurfaceClass(cbit.vcell.geometry.SurfaceClass) CSGeometry(org.sbml.jsbml.ext.spatial.CSGeometry) VCImage(cbit.image.VCImage) StructureMapping(cbit.vcell.mapping.StructureMapping) GeometryThumbnailImageFactoryAWT(cbit.vcell.geometry.GeometryThumbnailImageFactoryAWT) FeatureMapping(cbit.vcell.mapping.FeatureMapping) Structure(cbit.vcell.model.Structure) ModelUnitSystem(cbit.vcell.model.ModelUnitSystem) ParameterType(org.sbml.jsbml.ext.spatial.ParameterType) BioEventParameterType(cbit.vcell.mapping.BioEvent.BioEventParameterType) SampledFieldGeometry(org.sbml.jsbml.ext.spatial.SampledFieldGeometry) Geometry(cbit.vcell.geometry.Geometry) SampledFieldGeometry(org.sbml.jsbml.ext.spatial.SampledFieldGeometry) AnalyticGeometry(org.sbml.jsbml.ext.spatial.AnalyticGeometry) ParametricGeometry(org.sbml.jsbml.ext.spatial.ParametricGeometry) CSGeometry(org.sbml.jsbml.ext.spatial.CSGeometry) StringTokenizer(java.util.StringTokenizer) Expression(cbit.vcell.parser.Expression) Model(cbit.vcell.model.Model) BioModel(cbit.vcell.biomodel.BioModel) InterpolationKind(org.sbml.jsbml.ext.spatial.InterpolationKind) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) Parameter(org.sbml.jsbml.Parameter) ModelParameter(cbit.vcell.model.Model.ModelParameter) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) LocalParameter(org.sbml.jsbml.LocalParameter) KineticsProxyParameter(cbit.vcell.model.Kinetics.KineticsProxyParameter) UnresolvedParameter(cbit.vcell.model.Kinetics.UnresolvedParameter) CompartmentMapping(org.sbml.jsbml.ext.spatial.CompartmentMapping) Compartment(org.sbml.jsbml.Compartment) SpatialParameterPlugin(org.sbml.jsbml.ext.spatial.SpatialParameterPlugin) AnalyticGeometry(org.sbml.jsbml.ext.spatial.AnalyticGeometry) ExpressionException(cbit.vcell.parser.ExpressionException) GeometrySpec(cbit.vcell.geometry.GeometrySpec) DomainType(org.sbml.jsbml.ext.spatial.DomainType) SampledVolume(org.sbml.jsbml.ext.spatial.SampledVolume) ListOf(org.sbml.jsbml.ListOf) SpatialCompartmentPlugin(org.sbml.jsbml.ext.spatial.SpatialCompartmentPlugin) VCImageCompressed(cbit.image.VCImageCompressed) SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion) VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion) GeometricRegion(cbit.vcell.geometry.surface.GeometricRegion) AnalyticVolume(org.sbml.jsbml.ext.spatial.AnalyticVolume) InteriorPoint(org.sbml.jsbml.ext.spatial.InteriorPoint) ParametricGeometry(org.sbml.jsbml.ext.spatial.ParametricGeometry) SampledField(org.sbml.jsbml.ext.spatial.SampledField) Domain(org.sbml.jsbml.ext.spatial.Domain) GeometryClass(cbit.vcell.geometry.GeometryClass) GeometrySurfaceDescription(cbit.vcell.geometry.surface.GeometrySurfaceDescription) Extent(org.vcell.util.Extent) ISize(org.vcell.util.ISize) RegionInfo(cbit.vcell.geometry.RegionImage.RegionInfo) Membrane(cbit.vcell.model.Membrane) CSGObject(cbit.vcell.geometry.CSGObject) ImageSubVolume(cbit.vcell.geometry.ImageSubVolume) SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion) VCImageUncompressed(cbit.image.VCImageUncompressed) InteriorPoint(org.sbml.jsbml.ext.spatial.InteriorPoint) XMLStreamException(javax.xml.stream.XMLStreamException) SbmlException(org.vcell.sbml.SbmlException) IOException(java.io.IOException) PropertyVetoException(java.beans.PropertyVetoException) SBMLException(org.sbml.jsbml.SBMLException) ModelPropertyVetoException(cbit.vcell.model.ModelPropertyVetoException) ExpressionException(cbit.vcell.parser.ExpressionException)

Example 2 with AnalyticSubVolume

use of cbit.vcell.geometry.AnalyticSubVolume in project vcell by virtualcell.

the class XmlReader method getAnalyticSubVolume.

/**
 * This method returns an AnalyticSubVolume object from a XML representation.
 * Creation date: (5/1/2001 5:26:17 PM)
 * @return cbit.vcell.geometry.AnalyticSubVolume
 * @param param org.jdom.Element
 */
private AnalyticSubVolume getAnalyticSubVolume(Element param) throws XmlParseException {
    // retrieve the attributes
    String name = param.getAttributeValue(XMLTags.NameAttrTag);
    int handle = Integer.parseInt(param.getAttributeValue(XMLTags.HandleAttrTag));
    // process the key
    KeyValue key = null;
    String temp = param.getAttributeValue(XMLTags.KeyValueAttrTag);
    if (temp != null && temp.length() > 0 && this.readKeysFlag) {
        key = new KeyValue(temp);
    }
    // Retrieve the expression
    temp = param.getChildText(XMLTags.AnalyticExpressionTag, vcNamespace);
    if (temp == null) {
        throw new XmlParseException("A Problem occured while retrieving the analytic expression of the AnalyticSubvolume " + name);
    }
    Expression newexpression = unMangleExpression(temp);
    // Create the AnalyticCompartment
    AnalyticSubVolume newsubvolume = null;
    try {
        newsubvolume = new AnalyticSubVolume(key, name, newexpression, handle);
    } catch (ExpressionException e) {
        e.printStackTrace();
        throw new XmlParseException("An ExpressionException occured when creating the new AnalyticSubvolume " + name, e);
    }
    return newsubvolume;
}
Also used : KeyValue(org.vcell.util.document.KeyValue) Expression(cbit.vcell.parser.Expression) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) ExpressionException(cbit.vcell.parser.ExpressionException)

Example 3 with AnalyticSubVolume

use of cbit.vcell.geometry.AnalyticSubVolume in project vcell by virtualcell.

the class MathModel_SBMLExporter method addGeometry.

private static void addGeometry(Model sbmlModel, MathModel vcMathModel) {
    SpatialModelPlugin mplugin = (SpatialModelPlugin) sbmlModel.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
    // Creates a geometry object via SpatialModelPlugin object.
    org.sbml.jsbml.ext.spatial.Geometry sbmlGeometry = mplugin.getGeometry();
    sbmlGeometry.setCoordinateSystem(GeometryKind.cartesian);
    Geometry vcGeometry = vcMathModel.getGeometry();
    // 
    // list of CoordinateComponents : 1 if geometry is 1-d, 2 if geometry is 2-d, 3 if geometry is 3-d
    // 
    int dimension = vcGeometry.getDimension();
    Extent vcExtent = vcGeometry.getExtent();
    Origin vcOrigin = vcGeometry.getOrigin();
    // add x coordinate component
    CoordinateComponent coordCompX = sbmlGeometry.createCoordinateComponent();
    coordCompX.setSpatialId("CoordCompX");
    coordCompX.setType(CoordinateKind.cartesianX);
    Boundary minX = coordCompX.getBoundaryMaximum();
    minX.setSpatialId("Xmin");
    minX.setValue(vcOrigin.getX());
    Boundary maxX = coordCompX.getBoundaryMaximum();
    maxX.setSpatialId("Xmax");
    maxX.setValue(vcOrigin.getX() + (vcExtent.getX()));
    Parameter parameterX = sbmlModel.createParameter();
    // note for exporting BioModels rather than MathModels, get ReservedSymbol from Model with Role of ReservedSymbolRole.X
    parameterX.setId(ReservedVariable.X.getName());
    SpatialSymbolReference coordXSpatialRef = new SpatialSymbolReference();
    coordXSpatialRef.setSpatialRef(coordCompX.getSpatialId());
    SpatialParameterPlugin parameterXSpatialPlugin = (SpatialParameterPlugin) parameterX.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
    parameterXSpatialPlugin.setParamType(coordXSpatialRef);
    // add y coordinate component
    if (dimension == 2 || dimension == 3) {
        CoordinateComponent coordCompY = sbmlGeometry.createCoordinateComponent();
        coordCompY.setSpatialId("CoordCompY");
        coordCompY.setType(CoordinateKind.cartesianY);
        Boundary minY = coordCompY.getBoundaryMinimum();
        minY.setId("Ymin");
        minY.setValue(vcOrigin.getY());
        Boundary maxY = coordCompY.getBoundaryMaximum();
        maxY.setId("Ymax");
        maxY.setValue(vcOrigin.getY() + (vcExtent.getY()));
        Parameter parameterY = sbmlModel.createParameter();
        // note for exporting BioModels rather than MathModels, get ReservedSymbol from Model with Role of ReservedSymbolRole.Y
        parameterY.setId(ReservedVariable.Y.getName());
        SpatialSymbolReference coordYSpatialRef = new SpatialSymbolReference();
        coordYSpatialRef.setSpatialRef(coordCompY.getSpatialId());
        SpatialParameterPlugin parameterYSpatialPlugin = (SpatialParameterPlugin) parameterY.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
        parameterYSpatialPlugin.setParamType(coordYSpatialRef);
    }
    // add z coordinate component
    if (dimension == 3) {
        CoordinateComponent coordCompZ = sbmlGeometry.createCoordinateComponent();
        coordCompZ.setSpatialId("CoordCompZ");
        coordCompZ.setType(CoordinateKind.cartesianZ);
        Boundary minZ = coordCompZ.getBoundaryMinimum();
        minZ.setId("Zmin");
        minZ.setValue(vcOrigin.getZ());
        Boundary maxZ = coordCompZ.getBoundaryMaximum();
        maxZ.setId("Zmax");
        maxZ.setValue(vcOrigin.getZ() + (vcExtent.getZ()));
        Parameter parameterZ = sbmlModel.createParameter();
        // note for exporting BioModels rather than MathModels, get ReservedSymbol from Model with Role of ReservedSymbolRole.Y
        parameterZ.setId(ReservedVariable.Z.getName());
        SpatialSymbolReference coordZSpatialRef = new SpatialSymbolReference();
        coordZSpatialRef.setSpatialRef(coordCompZ.getSpatialId());
        SpatialParameterPlugin parameterZSpatialPlugin = (SpatialParameterPlugin) parameterZ.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
        parameterZSpatialPlugin.setParamType(coordZSpatialRef);
    }
    // 
    // list of domain types : subvolumes and surface classes from VC
    // Also create compartments - one compartment for each geometryClass. set id and spatialDimension based on type of geometryClass.
    // 
    boolean bAnalyticGeom = false;
    boolean bImageGeom = false;
    GeometryClass[] vcGeomClasses = vcGeometry.getGeometryClasses();
    int numVCGeomClasses = vcGeomClasses.length;
    for (int i = 0; i < numVCGeomClasses; i++) {
        DomainType domainType = sbmlGeometry.createDomainType();
        domainType.setId(vcGeomClasses[i].getName());
        if (vcGeomClasses[i] instanceof SubVolume) {
            if (((SubVolume) vcGeomClasses[i]) instanceof AnalyticSubVolume) {
                bAnalyticGeom = true;
            } else if (((SubVolume) vcGeomClasses[i]) instanceof ImageSubVolume) {
                bImageGeom = true;
            }
            domainType.setSpatialDimensions(3);
        } else if (vcGeomClasses[i] instanceof SurfaceClass) {
            domainType.setSpatialDimensions(2);
        }
    }
    // 
    // list of domains, adjacent domains : from VC geometricRegions
    // 
    GeometrySurfaceDescription vcGSD = vcGeometry.getGeometrySurfaceDescription();
    if (vcGSD.getRegionImage() == null) {
        try {
            vcGSD.updateAll();
        } catch (Exception e) {
            e.printStackTrace(System.out);
            throw new RuntimeException("Unable to generate region images for geometry");
        }
    }
    GeometricRegion[] vcGeometricRegions = vcGSD.getGeometricRegions();
    ISize sampleSize = vcGSD.getVolumeSampleSize();
    int numX = sampleSize.getX();
    int numY = sampleSize.getY();
    int numZ = sampleSize.getZ();
    double ox = vcOrigin.getX();
    double oy = vcOrigin.getY();
    double oz = vcOrigin.getZ();
    RegionInfo[] regionInfos = vcGSD.getRegionImage().getRegionInfos();
    Compartment compartment = null;
    for (int i = 0; i < vcGeometricRegions.length; i++) {
        // domains
        Domain domain = sbmlGeometry.createDomain();
        domain.setId(vcGeometricRegions[i].getName());
        compartment = sbmlModel.createCompartment();
        compartment.setId("compartment" + i);
        if (vcGeometricRegions[i] instanceof VolumeGeometricRegion) {
            domain.setDomainType(((VolumeGeometricRegion) vcGeometricRegions[i]).getSubVolume().getName());
            // domain.setImplicit(false);
            compartment.setSpatialDimensions(3);
            InteriorPoint interiorPt = domain.createInteriorPoint();
            int regionID = ((VolumeGeometricRegion) vcGeometricRegions[i]).getRegionID();
            boolean bFound = false;
            int regInfoIndx = 0;
            for (int j = 0; j < regionInfos.length; j++) {
                regInfoIndx = j;
                if (regionInfos[j].getRegionIndex() == regionID) {
                    int volIndx = 0;
                    for (int z = 0; z < numZ && !bFound; z++) {
                        for (int y = 0; y < numY && !bFound; y++) {
                            for (int x = 0; x < numX && !bFound; x++) {
                                if (regionInfos[j].isIndexInRegion(volIndx)) {
                                    bFound = true;
                                    double unit_z = (numZ > 1) ? ((double) z) / (numZ - 1) : 0.5;
                                    double coordZ = oz + vcExtent.getZ() * unit_z;
                                    double unit_y = (numY > 1) ? ((double) y) / (numY - 1) : 0.5;
                                    double coordY = oy + vcExtent.getY() * unit_y;
                                    double unit_x = (numX > 1) ? ((double) x) / (numX - 1) : 0.5;
                                    double coordX = ox + vcExtent.getX() * unit_x;
                                    interiorPt.setCoord1(coordX);
                                    interiorPt.setCoord2(coordY);
                                    interiorPt.setCoord3(coordZ);
                                }
                                volIndx++;
                            }
                        // end - for x
                        }
                    // end - for y
                    }
                // end - for z
                }
            // end if
            }
            // end for regionInfos
            if (!bFound) {
                throw new RuntimeException("Unable to find interior point for region '" + regionInfos[regInfoIndx].toString());
            }
        } else if (vcGeometricRegions[i] instanceof SurfaceGeometricRegion) {
            SurfaceGeometricRegion vcSurfaceGeomReg = (SurfaceGeometricRegion) vcGeometricRegions[i];
            GeometricRegion geomRegion0 = vcSurfaceGeomReg.getAdjacentGeometricRegions()[0];
            GeometricRegion geomRegion1 = vcSurfaceGeomReg.getAdjacentGeometricRegions()[1];
            SurfaceClass surfaceClass = vcGSD.getSurfaceClass(((VolumeGeometricRegion) geomRegion0).getSubVolume(), ((VolumeGeometricRegion) geomRegion1).getSubVolume());
            domain.setDomainType(surfaceClass.getName());
            // domain.setImplicit(true);
            compartment.setSpatialDimensions(2);
            // adjacent domains : 2 adjacent domain objects for each surfaceClass in VC.
            // adjacent domain 1
            AdjacentDomains adjDomain = sbmlGeometry.createAdjacentDomain();
            adjDomain.setId(TokenMangler.mangleToSName(vcSurfaceGeomReg.getName() + "_" + geomRegion0.getName()));
            adjDomain.setDomain1(vcSurfaceGeomReg.getName());
            adjDomain.setDomain2(geomRegion0.getName());
            // adjacent domain 2
            adjDomain = sbmlGeometry.createAdjacentDomain();
            adjDomain.setId(TokenMangler.mangleToSName(vcSurfaceGeomReg.getName() + "_" + geomRegion1.getName()));
            adjDomain.setDomain1(vcSurfaceGeomReg.getName());
            adjDomain.setDomain2(geomRegion1.getName());
        }
        // 
        // Mathmodel does not have structureMapping, hence creating compartmentMapping while creating domains.
        // @TODO : how to assign unitSize for compartmentMapping?
        // 
        SpatialCompartmentPlugin cplugin = (SpatialCompartmentPlugin) compartment.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
        CompartmentMapping compMapping = cplugin.getCompartmentMapping();
        String compMappingId = TokenMangler.mangleToSName(domain.getDomainType() + "_" + compartment.getId());
        compMapping.setId(compMappingId);
        compMapping.setDomainType(TokenMangler.mangleToSName(domain.getDomainType()));
    // try {
    // compMapping.setUnitSize(1.0);
    // } catch (ExpressionException e) {
    // e.printStackTrace(System.out);
    // throw new RuntimeException("Unable to create compartment mapping for structureMapping '" + compMapping.getId() +"' : " + e.getMessage());
    // }
    }
    AnalyticGeometry sbmlAnalyticGeom = null;
    SampledFieldGeometry sbmlSFGeom = null;
    // both image and analytic subvolumes?? == not handled in SBML at this time.
    if (bAnalyticGeom && !bImageGeom) {
        sbmlAnalyticGeom = sbmlGeometry.createAnalyticGeometry();
        sbmlAnalyticGeom.setId(TokenMangler.mangleToSName(vcGeometry.getName()));
    } else if (bImageGeom && !bAnalyticGeom) {
        // assuming image based geometry if not analytic geometry
        sbmlSFGeom = sbmlGeometry.createSampledFieldGeometry();
        sbmlSFGeom.setId(TokenMangler.mangleToSName(vcGeometry.getName()));
    } else if (bAnalyticGeom && bImageGeom) {
        throw new RuntimeException("Export to SBML of a combination of Image-based and Analytic geometries is not supported yet.");
    } else if (!bAnalyticGeom && !bImageGeom) {
        throw new RuntimeException("Unknown geometry type.");
    }
    // 
    for (int i = 0; i < vcGeomClasses.length; i++) {
        if (vcGeomClasses[i] instanceof AnalyticSubVolume) {
            // add analytiVols to sbmlAnalyticGeometry
            if (sbmlAnalyticGeom != null) {
                AnalyticVolume analyticVol = sbmlAnalyticGeom.createAnalyticVolume();
                analyticVol.setId(vcGeomClasses[i].getName());
                analyticVol.setDomainType(vcGeomClasses[i].getName());
                analyticVol.setFunctionType(FunctionKind.layered);
                analyticVol.setOrdinal(i);
                Expression expr = ((AnalyticSubVolume) vcGeomClasses[i]).getExpression();
                try {
                    String mathMLStr = ExpressionMathMLPrinter.getMathML(expr, true);
                    ASTNode mathMLNode = ASTNode.readMathMLFromString(mathMLStr);
                    analyticVol.setMath(mathMLNode);
                } catch (Exception e) {
                    e.printStackTrace(System.out);
                    throw new RuntimeException("Error converting VC subvolume expression to mathML" + e.getMessage());
                }
            } else {
                throw new RuntimeException("SBML AnalyticGeometry is null.");
            }
        } else if (vcGeomClasses[i] instanceof ImageSubVolume) {
            // add sampledVols to sbmlSFGeometry
            if (sbmlSFGeom != null) {
                SampledVolume sampledVol = sbmlSFGeom.createSampledVolume();
                sampledVol.setId(vcGeomClasses[i].getName());
                sampledVol.setDomainType(vcGeomClasses[i].getName());
                sampledVol.setSampledValue(((ImageSubVolume) vcGeomClasses[i]).getPixelValue());
            } else {
                throw new RuntimeException("SBML SampledFieldGeometry is null.");
            }
        }
    }
    if (sbmlSFGeom != null) {
        // add sampledField to sampledFieldGeometry
        SampledField sampledField = sbmlGeometry.createSampledField();
        VCImage vcImage = vcGeometry.getGeometrySpec().getImage();
        sampledField.setId(vcImage.getName());
        sampledField.setNumSamples1(vcImage.getNumX());
        if (vcImage.getNumY() > 1) {
            sampledField.setNumSamples2(vcImage.getNumY());
        }
        if (vcImage.getNumZ() > 1) {
            sampledField.setNumSamples3(vcImage.getNumZ());
        }
        sampledField.setInterpolationType(InterpolationKind.nearestneighbor);
        sampledField.setDataType(DataKind.UINT8);
        // add image from vcGeometrySpec to sampledField.
        try {
            StringBuffer sb = new StringBuffer();
            byte[] imagePixelsBytes = vcImage.getPixelsCompressed();
            for (int i = 0; i < imagePixelsBytes.length; i++) {
                int uint8_sample = ((int) imagePixelsBytes[i]) & 0xff;
                sb.append(uint8_sample + " ");
            }
            sampledField.setSamplesLength(vcImage.getNumXYZ());
            sampledField.setSamples(sb.toString().trim());
        } catch (ImageException e) {
            e.printStackTrace(System.out);
            throw new RuntimeException("Unable to export image from VCell to SBML : " + e.getMessage());
        }
    }
}
Also used : Origin(org.vcell.util.Origin) CompartmentMapping(org.sbml.jsbml.ext.spatial.CompartmentMapping) Compartment(org.sbml.jsbml.Compartment) SpatialParameterPlugin(org.sbml.jsbml.ext.spatial.SpatialParameterPlugin) AnalyticGeometry(org.sbml.jsbml.ext.spatial.AnalyticGeometry) Boundary(org.sbml.jsbml.ext.spatial.Boundary) DomainType(org.sbml.jsbml.ext.spatial.DomainType) SampledVolume(org.sbml.jsbml.ext.spatial.SampledVolume) SubVolume(cbit.vcell.geometry.SubVolume) ImageSubVolume(cbit.vcell.geometry.ImageSubVolume) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) SpatialCompartmentPlugin(org.sbml.jsbml.ext.spatial.SpatialCompartmentPlugin) CoordinateComponent(org.sbml.jsbml.ext.spatial.CoordinateComponent) VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion) VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion) SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion) GeometricRegion(cbit.vcell.geometry.surface.GeometricRegion) SpatialSymbolReference(org.sbml.jsbml.ext.spatial.SpatialSymbolReference) AnalyticVolume(org.sbml.jsbml.ext.spatial.AnalyticVolume) InteriorPoint(org.sbml.jsbml.ext.spatial.InteriorPoint) SampledField(org.sbml.jsbml.ext.spatial.SampledField) Domain(org.sbml.jsbml.ext.spatial.Domain) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) GeometryClass(cbit.vcell.geometry.GeometryClass) ImageException(cbit.image.ImageException) GeometrySurfaceDescription(cbit.vcell.geometry.surface.GeometrySurfaceDescription) Extent(org.vcell.util.Extent) SurfaceClass(cbit.vcell.geometry.SurfaceClass) ISize(org.vcell.util.ISize) RegionInfo(cbit.vcell.geometry.RegionImage.RegionInfo) VCImage(cbit.image.VCImage) ASTNode(org.sbml.jsbml.ASTNode) ImageSubVolume(cbit.vcell.geometry.ImageSubVolume) SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion) SampledFieldGeometry(org.sbml.jsbml.ext.spatial.SampledFieldGeometry) SpatialModelPlugin(org.sbml.jsbml.ext.spatial.SpatialModelPlugin) InteriorPoint(org.sbml.jsbml.ext.spatial.InteriorPoint) XMLStreamException(javax.xml.stream.XMLStreamException) ImageException(cbit.image.ImageException) SBMLException(org.sbml.jsbml.SBMLException) AdjacentDomains(org.sbml.jsbml.ext.spatial.AdjacentDomains) Geometry(cbit.vcell.geometry.Geometry) AnalyticGeometry(org.sbml.jsbml.ext.spatial.AnalyticGeometry) SampledFieldGeometry(org.sbml.jsbml.ext.spatial.SampledFieldGeometry) Expression(cbit.vcell.parser.Expression) Parameter(org.sbml.jsbml.Parameter)

Example 4 with AnalyticSubVolume

use of cbit.vcell.geometry.AnalyticSubVolume in project vcell by virtualcell.

the class GeometrySurfaceDescription method propertyChange.

/**
 * This method gets called when a bound property is changed.
 * @param evt A PropertyChangeEvent object describing the event source
 *   	and the property that has changed.
 */
public void propertyChange(java.beans.PropertyChangeEvent evt) {
    if (evt.getSource() == this && evt.getPropertyName().equals("volumeSampleSize")) {
        ISize oldValue = (ISize) evt.getOldValue();
        ISize newValue = (ISize) evt.getNewValue();
        if (!oldValue.compareEqual(newValue)) {
            try {
                // nobody listens to this, updateAll() will propagate changes
                getRegionImage0().setDirty();
                getSurfaceCollection0().setDirty();
                fieldGeometricRegions.setDirty();
            } catch (Exception e) {
                e.printStackTrace(System.out);
            }
        }
    }
    if (evt.getSource() == this && evt.getPropertyName().equals("filterCutoffFrequency")) {
        Double oldValue = (Double) evt.getOldValue();
        Double newValue = (Double) evt.getNewValue();
        if (!oldValue.equals(newValue)) {
            try {
                getSurfaceCollection0().setDirty();
                fieldGeometricRegions.setDirty();
            } catch (Exception e) {
                e.printStackTrace(System.out);
            }
        }
    }
    if (evt.getSource() == this && evt.getPropertyName().equals(PROPERTY_NAME_GEOMETRIC_REGIONS)) {
        refreshSurfaceClasses();
    }
    if (evt.getSource() == getGeometry().getGeometrySpec() && (evt.getPropertyName().equals("extent") || evt.getPropertyName().equals("origin"))) {
        Matchable oldExtentOrOrigin = (Matchable) evt.getOldValue();
        Matchable newExtentOrOrigin = (Matchable) evt.getNewValue();
        if (!Compare.isEqual(oldExtentOrOrigin, newExtentOrOrigin)) {
            try {
                // nobody listens to this, updateAll() will propagate changes
                getRegionImage0().setDirty();
                getSurfaceCollection0().setDirty();
                fieldGeometricRegions.setDirty();
            } catch (Exception e) {
                e.printStackTrace(System.out);
            }
        }
    }
    if (evt.getSource() instanceof AnalyticSubVolume && evt.getPropertyName().equals("expression")) {
        Expression oldExpression = (Expression) evt.getOldValue();
        Expression newExpression = (Expression) evt.getNewValue();
        if (!Compare.isEqual(oldExpression, newExpression)) {
            try {
                // nobody listens to this, updateAll() will propagate changes
                getRegionImage0().setDirty();
                getSurfaceCollection0().setDirty();
                fieldGeometricRegions.setDirty();
            } catch (Exception e) {
                e.printStackTrace(System.out);
            }
        }
    }
    if (evt.getSource() instanceof CSGObject && evt.getPropertyName().equals(CSGObject.PROPERTY_NAME_ROOT)) {
        try {
            // nobody listens to this, updateAll() will propagate changes
            getRegionImage0().setDirty();
            getSurfaceCollection0().setDirty();
            fieldGeometricRegions.setDirty();
        } catch (Exception e) {
            e.printStackTrace(System.out);
        }
    }
    if (evt.getSource() instanceof SubVolume && evt.getPropertyName().equals("name")) {
        String oldName = (String) evt.getOldValue();
        String newName = (String) evt.getNewValue();
        if (!Compare.isEqual(oldName, newName)) {
            try {
                fieldGeometricRegions.setDirty();
            } catch (Exception e) {
                e.printStackTrace(System.out);
            }
        }
    }
    if (evt.getSource() == getGeometry().getGeometrySpec() && evt.getPropertyName().equals("subVolumes")) {
        SubVolume[] oldValue = (SubVolume[]) evt.getOldValue();
        SubVolume[] newValue = (SubVolume[]) evt.getNewValue();
        // 
        for (int i = 0; oldValue != null && i < oldValue.length; i++) {
            oldValue[i].removePropertyChangeListener(this);
        }
        for (int i = 0; newValue != null && i < newValue.length; i++) {
            newValue[i].addPropertyChangeListener(this);
        }
        // 
        if (oldValue == null || newValue == null || !Compare.isEqualStrict(oldValue, newValue)) {
            try {
                // nobody listens to this, updateAll() will propagate changes
                getRegionImage0().setDirty();
                getSurfaceCollection0().setDirty();
                fieldGeometricRegions.setDirty();
            } catch (Exception e) {
                e.printStackTrace(System.out);
            }
        } else if (fieldGeometricRegions.getCurrentValue() != null && oldValue != newValue) {
            // 
            for (int i = 0; i < newValue.length; i++) {
                SubVolume subVolume = newValue[i];
                for (int j = 0; j < fieldGeometricRegions.getCurrentValue().length; j++) {
                    if (fieldGeometricRegions.getCurrentValue()[j] instanceof VolumeGeometricRegion) {
                        VolumeGeometricRegion volumeRegion = (VolumeGeometricRegion) fieldGeometricRegions.getCurrentValue()[j];
                        // 
                        if (volumeRegion.getSubVolume().compareEqual(subVolume) && volumeRegion.getSubVolume() != subVolume) {
                            volumeRegion.setSubVolume(subVolume);
                        }
                    }
                }
            }
        }
    }
}
Also used : Expression(cbit.vcell.parser.Expression) ISize(org.vcell.util.ISize) SubVolume(cbit.vcell.geometry.SubVolume) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) CSGObject(cbit.vcell.geometry.CSGObject) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) ImageException(cbit.image.ImageException) ExpressionException(cbit.vcell.parser.ExpressionException) PropertyVetoException(java.beans.PropertyVetoException) GeometryException(cbit.vcell.geometry.GeometryException) Matchable(org.vcell.util.Matchable)

Example 5 with AnalyticSubVolume

use of cbit.vcell.geometry.AnalyticSubVolume in project vcell by virtualcell.

the class MathTestingUtilities method constructExactMath.

/**
 * constructExactMath()
 *
 * take an equation of the form:
 *
 *     d A
 *     --- = F(A,t)
 *     d t
 *
 * and create a new equation with a known exact solution 'A_exact' by adding a forcing function R(t)
 *
 *     d A
 *     --- = F(A,t) + R(t)
 *     d t
 *
 * where:
 *
 *             d A_exact
 *      R(t) = --------- - F(A_exact,t)
 *                d t
 *
 * solving for R(t) is done analytically.
 *
 * Creation date: (1/21/2003 10:47:54 AM)
 * @return cbit.vcell.math.MathDescription
 * @param mathDesc cbit.vcell.math.MathDescription
 */
public static MathDescription constructExactMath(MathDescription mathDesc, java.util.Random random, ConstructedSolutionTemplate constructedSolutionTemplate) throws ExpressionException, MathException, MappingException {
    if (mathDesc.hasFastSystems()) {
        throw new RuntimeException("SolverTest.constructExactMath() suppport for fastSystems not yet implemented.");
    }
    MathDescription exactMath = null;
    try {
        exactMath = (MathDescription) BeanUtils.cloneSerializable(mathDesc);
        exactMath.setDescription("constructed exact solution from MathDescription (" + mathDesc.getName() + ")");
        exactMath.setName("exact from " + mathDesc.getName());
    } catch (Throwable e) {
        e.printStackTrace(System.out);
        throw new RuntimeException("error cloning MathDescription: " + e.getMessage());
    }
    // 
    // preload the VariableHash with existing Variables (and Constants,Functions,etc) and then sort all at once.
    // 
    VariableHash varHash = new VariableHash();
    Enumeration<Variable> enumVar = exactMath.getVariables();
    while (enumVar.hasMoreElements()) {
        varHash.addVariable(enumVar.nextElement());
    }
    java.util.Enumeration<SubDomain> subDomainEnum = exactMath.getSubDomains();
    while (subDomainEnum.hasMoreElements()) {
        SubDomain subDomain = subDomainEnum.nextElement();
        Domain domain = new Domain(subDomain);
        java.util.Enumeration<Equation> equationEnum = subDomain.getEquations();
        if (subDomain instanceof MembraneSubDomain) {
            MembraneSubDomain memSubDomain = (MembraneSubDomain) subDomain;
            AnalyticSubVolume insideAnalyticSubVolume = (AnalyticSubVolume) exactMath.getGeometry().getGeometrySpec().getSubVolume(memSubDomain.getInsideCompartment().getName());
            Function[] outwardNormalFunctions = getOutwardNormal(insideAnalyticSubVolume.getExpression(), "_" + insideAnalyticSubVolume.getName());
            for (int i = 0; i < outwardNormalFunctions.length; i++) {
                varHash.addVariable(outwardNormalFunctions[i]);
            }
        }
        while (equationEnum.hasMoreElements()) {
            Equation equation = equationEnum.nextElement();
            if (equation.getExactSolution() != null) {
                throw new RuntimeException("exact solution already exists");
            }
            Enumeration<Constant> origMathConstants = mathDesc.getConstants();
            if (equation instanceof OdeEquation) {
                OdeEquation odeEquation = (OdeEquation) equation;
                Expression substitutedRateExp = substituteWithExactSolution(odeEquation.getRateExpression(), (CompartmentSubDomain) subDomain, exactMath);
                SolutionTemplate solutionTemplate = constructedSolutionTemplate.getSolutionTemplate(equation.getVariable().getName(), subDomain.getName());
                String varName = odeEquation.getVariable().getName();
                String initName = null;
                while (origMathConstants.hasMoreElements()) {
                    Constant constant = origMathConstants.nextElement();
                    if (constant.getName().startsWith(varName + "_" + subDomain.getName() + DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONC_UNIT_PREFIX)) {
                        initName = constant.getName();
                    }
                }
                String exactName = varName + "_" + subDomain.getName() + "_exact";
                String errorName = varName + "_" + subDomain.getName() + "_error";
                String origRateName = "_" + varName + "_" + subDomain.getName() + "_origRate";
                String substitutedRateName = "_" + varName + "_" + subDomain.getName() + "_substitutedRate";
                String exactTimeDerivativeName = "_" + varName + "_" + subDomain.getName() + "_exact_dt";
                Expression exactExp = solutionTemplate.getTemplateExpression();
                Expression errorExp = new Expression(exactName + " - " + varName);
                Expression origRateExp = new Expression(odeEquation.getRateExpression());
                Expression exactTimeDerivativeExp = exactExp.differentiate("t").flatten();
                Expression newRate = new Expression(origRateName + " - " + substitutedRateName + " + " + exactTimeDerivativeName);
                Constant[] constants = solutionTemplate.getConstants();
                for (int i = 0; i < constants.length; i++) {
                    varHash.addVariable(constants[i]);
                }
                Expression initExp = new Expression(exactExp);
                initExp.substituteInPlace(new Expression("t"), new Expression(0.0));
                varHash.addVariable(new Function(initName, initExp.flatten(), domain));
                varHash.addVariable(new Function(exactName, exactExp, domain));
                varHash.addVariable(new Function(errorName, errorExp, domain));
                varHash.addVariable(new Function(exactTimeDerivativeName, exactTimeDerivativeExp, domain));
                varHash.addVariable(new Function(origRateName, origRateExp, domain));
                varHash.addVariable(new Function(substitutedRateName, substitutedRateExp, domain));
                odeEquation.setRateExpression(newRate);
                odeEquation.setInitialExpression(new Expression(initName));
                odeEquation.setExactSolution(new Expression(exactName));
            } else if (equation instanceof PdeEquation) {
                PdeEquation pdeEquation = (PdeEquation) equation;
                Expression substitutedRateExp = substituteWithExactSolution(pdeEquation.getRateExpression(), (CompartmentSubDomain) subDomain, exactMath);
                SolutionTemplate solutionTemplate = constructedSolutionTemplate.getSolutionTemplate(equation.getVariable().getName(), subDomain.getName());
                String varName = pdeEquation.getVariable().getName();
                String initName = null;
                while (origMathConstants.hasMoreElements()) {
                    Constant constant = origMathConstants.nextElement();
                    if (constant.getName().startsWith(varName + "_" + subDomain.getName() + DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONC_UNIT_PREFIX)) {
                        initName = constant.getName();
                    }
                }
                String diffusionRateName = "_" + varName + "_" + subDomain.getName() + "_diffusionRate";
                String exactName = varName + "_" + subDomain.getName() + "_exact";
                String errorName = varName + "_" + subDomain.getName() + "_error";
                String origRateName = "_" + varName + "_" + subDomain.getName() + "_origRate";
                String substitutedRateName = "_" + varName + "_" + subDomain.getName() + "_substitutedRate";
                String exactTimeDerivativeName = "_" + varName + "_" + subDomain.getName() + "_exact_dt";
                String exactDxName = "_" + varName + "_" + subDomain.getName() + "_exact_dx";
                String exactDyName = "_" + varName + "_" + subDomain.getName() + "_exact_dy";
                String exactDzName = "_" + varName + "_" + subDomain.getName() + "_exact_dz";
                String exactDx2Name = "_" + varName + "_" + subDomain.getName() + "_exact_dx2";
                String exactDy2Name = "_" + varName + "_" + subDomain.getName() + "_exact_dy2";
                String exactDz2Name = "_" + varName + "_" + subDomain.getName() + "_exact_dz2";
                String exactLaplacianName = "_" + varName + "_" + subDomain.getName() + "_exact_laplacian";
                Expression exactExp = solutionTemplate.getTemplateExpression();
                Expression errorExp = new Expression(exactName + " - " + varName);
                Expression origRateExp = new Expression(pdeEquation.getRateExpression());
                Expression initExp = new Expression(exactExp);
                initExp.substituteInPlace(new Expression("t"), new Expression(0.0));
                initExp = initExp.flatten();
                Expression exactTimeDerivativeExp = exactExp.differentiate("t").flatten();
                Expression exactDxExp = exactExp.differentiate("x").flatten();
                Expression exactDx2Exp = exactDxExp.differentiate("x").flatten();
                Expression exactDyExp = exactExp.differentiate("y").flatten();
                Expression exactDy2Exp = exactDxExp.differentiate("y").flatten();
                Expression exactDzExp = exactExp.differentiate("z").flatten();
                Expression exactDz2Exp = exactDxExp.differentiate("z").flatten();
                Expression exactLaplacianExp = Expression.add(Expression.add(exactDx2Exp, exactDy2Exp), exactDz2Exp).flatten();
                Expression newRate = new Expression(origRateName + " - " + substitutedRateName + " - ((" + diffusionRateName + ")*" + exactLaplacianName + ")" + " + " + exactTimeDerivativeName);
                Constant[] constants = solutionTemplate.getConstants();
                for (int i = 0; i < constants.length; i++) {
                    varHash.addVariable(constants[i]);
                }
                varHash.addVariable(new Function(initName, initExp, domain));
                varHash.addVariable(new Function(diffusionRateName, new Expression(pdeEquation.getDiffusionExpression()), domain));
                varHash.addVariable(new Function(exactName, exactExp, domain));
                varHash.addVariable(new Function(errorName, errorExp, domain));
                varHash.addVariable(new Function(exactTimeDerivativeName, exactTimeDerivativeExp, domain));
                varHash.addVariable(new Function(origRateName, origRateExp, domain));
                varHash.addVariable(new Function(substitutedRateName, substitutedRateExp, domain));
                varHash.addVariable(new Function(exactDxName, exactDxExp, domain));
                varHash.addVariable(new Function(exactDyName, exactDyExp, domain));
                varHash.addVariable(new Function(exactDzName, exactDzExp, domain));
                varHash.addVariable(new Function(exactDx2Name, exactDx2Exp, domain));
                varHash.addVariable(new Function(exactDy2Name, exactDy2Exp, domain));
                varHash.addVariable(new Function(exactDz2Name, exactDz2Exp, domain));
                varHash.addVariable(new Function(exactLaplacianName, exactLaplacianExp, domain));
                pdeEquation.setRateExpression(newRate);
                pdeEquation.setInitialExpression(new Expression(initName));
                pdeEquation.setDiffusionExpression(new Expression(diffusionRateName));
                pdeEquation.setExactSolution(new Expression(exactName));
                CompartmentSubDomain compartmentSubDomain = (CompartmentSubDomain) subDomain;
                if (compartmentSubDomain.getBoundaryConditionXm().isDIRICHLET()) {
                    Expression origExp = pdeEquation.getBoundaryXm();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryXm(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
                    } else {
                        pdeEquation.setBoundaryXm(exactExp);
                    }
                } else if (compartmentSubDomain.getBoundaryConditionXm().isNEUMANN()) {
                    Expression origExp = pdeEquation.getBoundaryXm();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryXm(new Expression(origExp + "-" + substitutedExp + "-" + diffusionRateName + "*" + exactDxName));
                    } else {
                        pdeEquation.setBoundaryXm(new Expression("-" + diffusionRateName + "*" + exactDxName));
                    }
                } else {
                    throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionXm());
                }
                if (compartmentSubDomain.getBoundaryConditionXp().isDIRICHLET()) {
                    Expression origExp = pdeEquation.getBoundaryXp();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryXp(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
                    } else {
                        pdeEquation.setBoundaryXp(exactExp);
                    }
                } else if (compartmentSubDomain.getBoundaryConditionXp().isNEUMANN()) {
                    Expression origExp = pdeEquation.getBoundaryXp();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryXp(new Expression(origExp + "-" + substitutedExp + "+" + diffusionRateName + "*" + exactDxName));
                    } else {
                        pdeEquation.setBoundaryXp(new Expression(diffusionRateName + "*" + exactDxName));
                    }
                } else {
                    throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionXp());
                }
                if (compartmentSubDomain.getBoundaryConditionYm().isDIRICHLET()) {
                    Expression origExp = pdeEquation.getBoundaryYm();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryYm(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
                    } else {
                        pdeEquation.setBoundaryYm(exactExp);
                    }
                } else if (compartmentSubDomain.getBoundaryConditionYm().isNEUMANN()) {
                    Expression origExp = pdeEquation.getBoundaryYm();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryYm(new Expression(origExp + "-" + substitutedExp + "-" + diffusionRateName + "*" + exactDyName));
                    } else {
                        pdeEquation.setBoundaryYm(new Expression("-" + diffusionRateName + "*" + exactDyName));
                    }
                } else {
                    throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionYm());
                }
                if (compartmentSubDomain.getBoundaryConditionYp().isDIRICHLET()) {
                    Expression origExp = pdeEquation.getBoundaryYp();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryYp(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
                    } else {
                        pdeEquation.setBoundaryYp(exactExp);
                    }
                } else if (compartmentSubDomain.getBoundaryConditionYp().isNEUMANN()) {
                    Expression origExp = pdeEquation.getBoundaryYp();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryYp(new Expression(origExp + "-" + substitutedExp + "+" + diffusionRateName + "*" + exactDyName));
                    } else {
                        pdeEquation.setBoundaryYp(new Expression(diffusionRateName + "*" + exactDyName));
                    }
                } else {
                    throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionYp());
                }
                if (compartmentSubDomain.getBoundaryConditionZm().isDIRICHLET()) {
                    Expression origExp = pdeEquation.getBoundaryZm();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryZm(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
                    } else {
                        pdeEquation.setBoundaryZm(exactExp);
                    }
                } else if (compartmentSubDomain.getBoundaryConditionZm().isNEUMANN()) {
                    Expression origExp = pdeEquation.getBoundaryZm();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryZm(new Expression(origExp + "-" + substitutedExp + "-" + diffusionRateName + "*" + exactDzName));
                    } else {
                        pdeEquation.setBoundaryZm(new Expression("-" + diffusionRateName + "*" + exactDzName));
                    }
                } else {
                    throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionXm());
                }
                if (compartmentSubDomain.getBoundaryConditionZp().isDIRICHLET()) {
                    Expression origExp = pdeEquation.getBoundaryZp();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryZp(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
                    } else {
                        pdeEquation.setBoundaryZp(exactExp);
                    }
                } else if (compartmentSubDomain.getBoundaryConditionZp().isNEUMANN()) {
                    Expression origExp = pdeEquation.getBoundaryZp();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryZp(new Expression(origExp + "-" + substitutedExp + "+" + diffusionRateName + "*" + exactDzName));
                    } else {
                        pdeEquation.setBoundaryZp(new Expression(diffusionRateName + "*" + exactDzName));
                    }
                } else {
                    throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionZp());
                }
            } else {
                throw new RuntimeException("SolverTest.constructedExactMath(): equation type " + equation.getClass().getName() + " not yet implemented");
            }
        }
        if (subDomain instanceof MembraneSubDomain) {
            MembraneSubDomain membraneSubDomain = (MembraneSubDomain) subDomain;
            Enumeration<JumpCondition> enumJumpConditions = membraneSubDomain.getJumpConditions();
            while (enumJumpConditions.hasMoreElements()) {
                JumpCondition jumpCondition = enumJumpConditions.nextElement();
                Expression origInfluxExp = jumpCondition.getInFluxExpression();
                Expression origOutfluxExp = jumpCondition.getOutFluxExpression();
                Expression substitutedInfluxExp = substituteWithExactSolution(origInfluxExp, membraneSubDomain, exactMath);
                Expression substitutedOutfluxExp = substituteWithExactSolution(origOutfluxExp, membraneSubDomain, exactMath);
                String varName = jumpCondition.getVariable().getName();
                String origInfluxName = "_" + varName + "_" + subDomain.getName() + "_origInflux";
                String origOutfluxName = "_" + varName + "_" + subDomain.getName() + "_origOutflux";
                String substitutedInfluxName = "_" + varName + "_" + subDomain.getName() + "_substitutedInflux";
                String substitutedOutfluxName = "_" + varName + "_" + subDomain.getName() + "_substitutedOutflux";
                String diffusionRateInsideName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_diffusionRate";
                String diffusionRateOutsideName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_diffusionRate";
                String exactInsideDxName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_exact_dx";
                String exactInsideDyName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_exact_dy";
                String exactInsideDzName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_exact_dz";
                String exactOutsideDxName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_exact_dx";
                String exactOutsideDyName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_exact_dy";
                String exactOutsideDzName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_exact_dz";
                String outwardNormalXName = "_" + membraneSubDomain.getInsideCompartment().getName() + "_Nx";
                String outwardNormalYName = "_" + membraneSubDomain.getInsideCompartment().getName() + "_Ny";
                String outwardNormalZName = "_" + membraneSubDomain.getInsideCompartment().getName() + "_Nz";
                String exactInfluxName = "_" + varName + "_" + membraneSubDomain.getName() + "_exactInflux";
                String exactOutfluxName = "_" + varName + "_" + membraneSubDomain.getName() + "_exactOutflux";
                Expression exactInfluxExp = new Expression(diffusionRateInsideName + " * (" + outwardNormalXName + "*" + exactInsideDxName + " + " + outwardNormalYName + "*" + exactInsideDyName + " + " + outwardNormalZName + "*" + exactInsideDzName + ")");
                Expression exactOutfluxExp = new Expression("-" + diffusionRateOutsideName + " * (" + outwardNormalXName + "*" + exactOutsideDxName + " + " + outwardNormalYName + "*" + exactOutsideDyName + " + " + outwardNormalZName + "*" + exactOutsideDzName + ")");
                Expression newInfluxExp = new Expression(origInfluxName + " - " + substitutedInfluxName + " + " + exactInfluxName);
                Expression newOutfluxExp = new Expression(origOutfluxName + " - " + substitutedOutfluxName + " + " + exactOutfluxName);
                varHash.addVariable(new Function(origInfluxName, origInfluxExp, domain));
                varHash.addVariable(new Function(origOutfluxName, origOutfluxExp, domain));
                varHash.addVariable(new Function(exactInfluxName, exactInfluxExp, domain));
                varHash.addVariable(new Function(exactOutfluxName, exactOutfluxExp, domain));
                varHash.addVariable(new Function(substitutedInfluxName, substitutedInfluxExp, domain));
                varHash.addVariable(new Function(substitutedOutfluxName, substitutedOutfluxExp, domain));
                jumpCondition.setInFlux(newInfluxExp);
                jumpCondition.setOutFlux(newOutfluxExp);
            }
        }
    }
    exactMath.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
    if (!exactMath.isValid()) {
        throw new RuntimeException("generated Math is not valid: " + exactMath.getWarning());
    }
    return exactMath;
}
Also used : JumpCondition(cbit.vcell.math.JumpCondition) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable) InsideVariable(cbit.vcell.math.InsideVariable) SensVariable(cbit.vcell.solver.ode.SensVariable) FilamentVariable(cbit.vcell.math.FilamentVariable) VolVariable(cbit.vcell.math.VolVariable) VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable) ReservedVariable(cbit.vcell.math.ReservedVariable) MemVariable(cbit.vcell.math.MemVariable) OutsideVariable(cbit.vcell.math.OutsideVariable) FilamentRegionVariable(cbit.vcell.math.FilamentRegionVariable) Variable(cbit.vcell.math.Variable) MathDescription(cbit.vcell.math.MathDescription) VariableHash(cbit.vcell.math.VariableHash) Constant(cbit.vcell.math.Constant) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) PdeEquation(cbit.vcell.math.PdeEquation) Function(cbit.vcell.math.Function) OdeEquation(cbit.vcell.math.OdeEquation) PdeEquation(cbit.vcell.math.PdeEquation) Equation(cbit.vcell.math.Equation) SolutionTemplate(cbit.vcell.numericstest.SolutionTemplate) ConstructedSolutionTemplate(cbit.vcell.numericstest.ConstructedSolutionTemplate) OdeEquation(cbit.vcell.math.OdeEquation) Expression(cbit.vcell.parser.Expression) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) Domain(cbit.vcell.math.Variable.Domain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume)

Aggregations

AnalyticSubVolume (cbit.vcell.geometry.AnalyticSubVolume)20 Expression (cbit.vcell.parser.Expression)14 CSGObject (cbit.vcell.geometry.CSGObject)13 SubVolume (cbit.vcell.geometry.SubVolume)11 GeometryThumbnailImageFactoryAWT (cbit.vcell.geometry.GeometryThumbnailImageFactoryAWT)10 Geometry (cbit.vcell.geometry.Geometry)8 ImageSubVolume (cbit.vcell.geometry.ImageSubVolume)7 Extent (org.vcell.util.Extent)7 AsynchClientTask (cbit.vcell.client.task.AsynchClientTask)6 Hashtable (java.util.Hashtable)6 Origin (org.vcell.util.Origin)6 GeometrySpec (cbit.vcell.geometry.GeometrySpec)5 ExpressionException (cbit.vcell.parser.ExpressionException)5 ISize (org.vcell.util.ISize)5 ImageException (cbit.image.ImageException)4 VCImage (cbit.image.VCImage)4 BioModel (cbit.vcell.biomodel.BioModel)4 GeometrySurfaceDescription (cbit.vcell.geometry.surface.GeometrySurfaceDescription)4 SimulationContext (cbit.vcell.mapping.SimulationContext)4 GeometryClass (cbit.vcell.geometry.GeometryClass)3