Examples with CompartmentSubDomain cbit.vcell.math.CompartmentSubDomain used on opensource projects

Example 31 with CompartmentSubDomain

use of cbit.vcell.math.CompartmentSubDomain in project vcell by virtualcell.

the class XmlReader method getMembraneSubDomain.

/**
 * This method returns a MembraneSubDomain object from a XML Element.
 * Creation date: (5/18/2001 4:23:30 PM)
 * @return cbit.vcell.math.MembraneSubDomain
 * @param param org.jdom.Element
 * @exception cbit.vcell.xml.XmlParseException The exception description.
 */
@SuppressWarnings("unchecked")
private MembraneSubDomain getMembraneSubDomain(Element param, MathDescription mathDesc) throws XmlParseException {
    // no need to do anything with the 'Name' attribute : constructor of MembraneSubDomain creates name from inside/outside compartmentSubDomains.
    // String msdName = unMangle( param.getAttributeValue(XMLTags.NameAttrTag) );
    // if ( msdName != null) {
    // }
    // get compartmentSubDomain references
    // inside
    String name = unMangle(param.getAttributeValue(XMLTags.InsideCompartmentTag));
    CompartmentSubDomain insideRef = (CompartmentSubDomain) mathDesc.getCompartmentSubDomain(name);
    if (insideRef == null) {
        throw new XmlParseException("The reference to the inside CompartmentSubDomain " + name + ", could not be resolved!");
    }
    // outside
    name = unMangle(param.getAttributeValue(XMLTags.OutsideCompartmentTag));
    CompartmentSubDomain outsideRef = (CompartmentSubDomain) mathDesc.getCompartmentSubDomain(name);
    if (outsideRef == null) {
        throw new XmlParseException("The reference to the outside CompartmentSubDomain " + name + ", could not be resolved!");
    }
    // *** create new Membrane SubDomain ***
    SubVolume insideSubVolume = mathDesc.getGeometry().getGeometrySpec().getSubVolume(insideRef.getName());
    SubVolume outsideSubVolume = mathDesc.getGeometry().getGeometrySpec().getSubVolume(outsideRef.getName());
    SurfaceClass surfaceClass = mathDesc.getGeometry().getGeometrySurfaceDescription().getSurfaceClass(insideSubVolume, outsideSubVolume);
    MembraneSubDomain subDomain = new MembraneSubDomain(insideRef, outsideRef, surfaceClass.getName());
    transcribeComments(param, subDomain);
    // Process BoundaryConditions
    Iterator<Element> iterator = param.getChildren(XMLTags.BoundaryTypeTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempelement = (Element) iterator.next();
        // create BoundaryConditionType
        String temp = tempelement.getAttributeValue(XMLTags.BoundaryTypeAttrTag);
        BoundaryConditionType bType = new BoundaryConditionType(temp);
        // Process Xm
        if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueXm)) {
            subDomain.setBoundaryConditionXm(bType);
        } else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueXp)) {
            // Process Xp
            subDomain.setBoundaryConditionXp(bType);
        } else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueYm)) {
            // Process Ym
            subDomain.setBoundaryConditionYm(bType);
        } else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueYp)) {
            // Process Yp
            subDomain.setBoundaryConditionYp(bType);
        } else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueZm)) {
            // Process Zm
            subDomain.setBoundaryConditionZm(bType);
        } else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueZp)) {
            // Process Zp
            subDomain.setBoundaryConditionZp(bType);
        } else {
            // If not indentified throw an exception!!
            throw new XmlParseException("Unknown BoundaryConditionType: " + tempelement.getAttributeValue(XMLTags.BoundaryAttrTag));
        }
    }
    // Add OdeEquations
    iterator = param.getChildren(XMLTags.OdeEquationTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempElement = (Element) iterator.next();
        OdeEquation odeEquation = getOdeEquation(tempElement, mathDesc);
        try {
            subDomain.addEquation(odeEquation);
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding an OdeEquation to a MembraneSubDomain!", e);
        }
    }
    // process PdeEquations
    iterator = param.getChildren(XMLTags.PdeEquationTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempElement = (Element) iterator.next();
        try {
            subDomain.addEquation(getPdeEquation(tempElement, mathDesc));
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding an PdeEquation to the MembraneSubDomain " + name, e);
        }
    }
    // Add JumpConditions
    iterator = param.getChildren(XMLTags.JumpConditionTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempElement = (Element) iterator.next();
        try {
            subDomain.addJumpCondition(getJumpCondition(tempElement, mathDesc));
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding a JumpCondition to a MembraneSubDomain!", e);
        }
    }
    // Add the FastSystem (if any)
    Element tempElement = param.getChild(XMLTags.FastSystemTag, vcNamespace);
    if (tempElement != null) {
        subDomain.setFastSystem(getFastSystem(tempElement, mathDesc));
    }
    // add MembraneRegionEquation
    iterator = param.getChildren(XMLTags.MembraneRegionEquationTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        tempElement = (Element) iterator.next();
        try {
            subDomain.addEquation(getMembraneRegionEquation(tempElement, mathDesc));
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding a MembraneRegionEquation to a MEmbraneSubDomain!", e);
        }
    }
    iterator = param.getChildren(XMLTags.ParticleJumpProcessTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempelement = (Element) iterator.next();
        try {
            subDomain.addParticleJumpProcess(getParticleJumpProcess(tempelement, mathDesc));
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding a jump process to the MembraneSubDomain " + name, e);
        }
    }
    iterator = param.getChildren(XMLTags.ParticlePropertiesTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempelement = (Element) iterator.next();
        try {
            subDomain.addParticleProperties(getParticleProperties(tempelement, mathDesc));
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding a jump process to the MembraneSubDomain " + name, e);
        }
    }
    // process ComputeNormal "equations"
    iterator = param.getChildren(XMLTags.ComputeNormalTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempelement = (Element) iterator.next();
        try {
            subDomain.addEquation(getComputeNormal(tempelement, mathDesc));
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding an ComputeNormal 'equation' to the MembraneSubDomain " + name, e);
        }
    }
    Element velElem = param.getChild(XMLTags.VelocityTag, vcNamespace);
    setMembraneSubdomainVelocity(velElem, XMLTags.XAttrTag, subDomain::setVelocityX);
    setMembraneSubdomainVelocity(velElem, XMLTags.YAttrTag, subDomain::setVelocityY);
    return subDomain;
}

Example 32 with CompartmentSubDomain

use of cbit.vcell.math.CompartmentSubDomain in project vcell by virtualcell.

the class MathTestingUtilities method comparePDEResults.

/**
 * Insert the method's description here.
 * Creation date: (8/20/2003 12:58:10 PM)
 */
public static SimulationComparisonSummary comparePDEResults(SimulationSymbolTable testSimSymbolTable, PDEDataManager testDataManager, SimulationSymbolTable refSimSymbolTable, PDEDataManager refDataManager, String[] varsToCompare, double absErrorThreshold, double relErrorThreshold, VCDocument refDocument, DataInfoProvider refDataInfoProvider, VCDocument testDocument, DataInfoProvider testDataInfoProvider) throws DataAccessException, ExpressionException {
    java.util.Hashtable<String, DataErrorSummary> tempVarHash = new java.util.Hashtable<String, DataErrorSummary>();
    boolean bTimesEqual = true;
    double[] testTimeArray = testDataManager.getDataSetTimes();
    double[] refTimeArray = refDataManager.getDataSetTimes();
    if (testTimeArray.length != refTimeArray.length) {
        bTimesEqual = false;
    // throw new RuntimeException("Data times for reference and test simulations don't match, cannot compare the two simulations!");
    } else {
        for (int i = 0; i < testTimeArray.length; i++) {
            if (testTimeArray[i] != refTimeArray[i]) {
                bTimesEqual = false;
            }
        }
    }
    // if (!checkSimVars(testSim, refSim)) {
    // //String errorS = "TestVars - ";
    // Variable[] testVars = testSim.getVariables();
    // //for(int i =0;i<vars.length;i+= 1){
    // //errorS+= vars[i].getName()+" ";
    // //}
    // //errorS+=" <<>> RefVars - ";
    // Variable[] refVars = refSim.getVariables();
    // //for(int i =0;i<vars.length;i+= 1){
    // //errorS+= vars[i].getName()+" ";
    // //}
    // throw new RuntimeException(
    // "VarNotMatch testLength="+(testVars != null?testVars.length+"":"null")+" refLength="+(refVars != null?refVars.length+"":"null"));
    // }
    CartesianMesh testMesh = testDataManager.getMesh();
    MathDescription testMathDesc = testSimSymbolTable.getSimulation().getMathDescription();
    // Variable[] refVars = refSim.getVariables();
    MathDescription refMathDesc = refSimSymbolTable.getSimulation().getMathDescription();
    CartesianMesh refMesh = refDataManager.getMesh();
    int[] membraneIndexMapping = null;
    // Get volumeSubdomains from mathDesc/mesh and store in lookupTable for testSimulation
    int testNumVol = testMesh.getSizeX() * testMesh.getSizeY() * testMesh.getSizeZ();
    CompartmentSubDomain[] testVolSubDomainLookup = new CompartmentSubDomain[testNumVol];
    for (int i = 0; i < testNumVol; i++) {
        int subVolumeIndex = testMesh.getSubVolumeFromVolumeIndex(i);
        SubVolume subVolume = testMathDesc.getGeometry().getGeometrySpec().getSubVolume(subVolumeIndex);
        CompartmentSubDomain compSubDomain = testMathDesc.getCompartmentSubDomain(subVolume.getName());
        testVolSubDomainLookup[i] = compSubDomain;
    }
    // Get membraneSubdomains from mathDesc/mesh and store in lookupTable for testSimulation
    int testNumMem = testMesh.getMembraneElements().length;
    MembraneSubDomain[] testMemSubDomainLookup = new MembraneSubDomain[testNumMem];
    for (int i = 0; i < testNumMem; i++) {
        int insideVolIndex = testMesh.getMembraneElements()[i].getInsideVolumeIndex();
        int outsideVolIndex = testMesh.getMembraneElements()[i].getOutsideVolumeIndex();
        MembraneSubDomain memSubDomain = testMathDesc.getMembraneSubDomain(testVolSubDomainLookup[insideVolIndex], testVolSubDomainLookup[outsideVolIndex]);
        testMemSubDomainLookup[i] = memSubDomain;
    }
    // Get volumeSubdomains from mathDesc/mesh and store in lookupTable for refSimulation
    int refNumVol = refMesh.getSizeX() * refMesh.getSizeY() * refMesh.getSizeZ();
    CompartmentSubDomain[] refVolSubDomainLookup = new CompartmentSubDomain[refNumVol];
    for (int i = 0; i < refNumVol; i++) {
        int subVolumeIndex = refMesh.getSubVolumeFromVolumeIndex(i);
        SubVolume subVolume = refMathDesc.getGeometry().getGeometrySpec().getSubVolume(subVolumeIndex);
        CompartmentSubDomain compSubDomain = refMathDesc.getCompartmentSubDomain(subVolume.getName());
        refVolSubDomainLookup[i] = compSubDomain;
    }
    // Get membraneSubdomains from mathDesc/mesh and store in lookupTable for refSimulation
    int refNumMem = refMesh.getMembraneElements().length;
    MembraneSubDomain[] refMemSubDomainLookup = new MembraneSubDomain[refNumMem];
    for (int i = 0; i < refNumMem; i++) {
        int insideVolIndex = refMesh.getMembraneElements()[i].getInsideVolumeIndex();
        int outsideVolIndex = refMesh.getMembraneElements()[i].getOutsideVolumeIndex();
        MembraneSubDomain memSubDomain = refMathDesc.getMembraneSubDomain(refVolSubDomainLookup[insideVolIndex], refVolSubDomainLookup[outsideVolIndex]);
        refMemSubDomainLookup[i] = memSubDomain;
    }
    SimulationComparisonSummary simComparisonSummary = new SimulationComparisonSummary();
    String hashKey = new String("");
    DataErrorSummary tempVar = null;
    DataIdentifier[] refDataIDs = refDataManager.getDataIdentifiers();
    // for each var, do the following :
    for (int i = 0; i < varsToCompare.length; i++) {
        DataIdentifier refDataID = null;
        for (int j = 0; j < refDataIDs.length; j++) {
            if (refDataIDs[j].getName().equals(varsToCompare[i])) {
                refDataID = refDataIDs[j];
                break;
            }
        }
        // 
        // Find REFERENCE variable
        // 
        Variable refVar = getSimVar(refSimSymbolTable, varsToCompare[i]);
        if (refVar == null) {
            // Should only happen if TEST sims were generated 'post-domains' and REFERENCE sims were generated 'pre-domains'
            if (refDataID != null) {
                throw new RuntimeException("Unexpected reference condition: '" + varsToCompare[i] + "' not found in symboltable but was found in dataidentifiers");
            }
            if (testDocument instanceof BioModel) {
                // Only BioModels need to be checked
                // Look in TEST for a speciescontext with matching species name
                System.out.println("ReferenceVariable: using alternate method to find '" + varsToCompare[i] + "'");
                BioModel refBioModel = (BioModel) refDocument;
                BioModel testBioModel = (BioModel) testDocument;
                SpeciesContext testSpeciesContext = testBioModel.getModel().getSpeciesContext(varsToCompare[i]);
                if (testSpeciesContext != null) {
                    refVar = refSimSymbolTable.getVariable(testSpeciesContext.getSpecies().getCommonName());
                    for (int j = 0; j < refDataIDs.length; j++) {
                        if (refDataIDs[j].getName().equals(testSpeciesContext.getSpecies().getCommonName())) {
                            refDataID = refDataIDs[j];
                            break;
                        }
                    }
                }
            }
            if (refVar == null || refDataID == null) {
                Simulation refSim = refSimSymbolTable.getSimulation();
                throw new RuntimeException("The variable " + varsToCompare[i] + " was not found in Simulation (" + refSim.getName() + " " + refSim.getVersion().getDate() + ")\n");
            }
        }
        // 
        // Find TEST variable (assumed to have sims generated with a software version later than REFERENCE)
        // 
        Variable testVar = getSimVar(testSimSymbolTable, varsToCompare[i]);
        if (testVar == null) {
            // Should only happen if TEST sims were generated 'post-domains' and REFERENCE sims were generated 'pre-domains'
            System.out.println("TestVariable: using alternate method to find '" + varsToCompare[i] + "'");
            BioModel testBioModel = (BioModel) testDocument;
            SpeciesContext[] speciesContexts = testBioModel.getModel().getSpeciesContexts();
            boolean bSkip = false;
            for (int j = 0; j < speciesContexts.length; j++) {
                if (speciesContexts[j].getSpecies().getCommonName().equals(varsToCompare[i])) {
                    testVar = testSimSymbolTable.getVariable(speciesContexts[j].getName());
                    if (testVar == null) {
                        throw new RuntimeException("Speciescontext name '" + speciesContexts[j].getName() + "' not found in testsimsymboltable");
                    }
                    // If we got here it means at least one matching speciescontext was found in TEST with
                    // a species name matching varsToCompare[i].  We can skip because the matching speciesconetext
                    // will be used to do a comparison at some point.
                    bSkip = true;
                    break;
                }
            }
            if (bSkip) {
                // these are tested already using full simcontext names
                System.out.println("Skipping '" + varsToCompare[i] + "' as lookup in testSimSymbolTable");
                continue;
            }
            Simulation refSim = refSimSymbolTable.getSimulation();
            throw new RuntimeException("The variable " + varsToCompare[i] + " was not found in Simulation (" + refSim.getName() + " " + refSim.getVersion().getDate() + ")\n");
        }
        // for each time in timeArray. ('t' is used to index the testTimeArray, for interpolation purposes.)
        int t = 0;
        for (int j = 0; j < refTimeArray.length; j++) {
            // get data block from varName, data from datablock
            SimDataBlock refSimDataBlock = refDataManager.getSimDataBlock(refVar.getName(), refTimeArray[j]);
            double[] refData = refSimDataBlock.getData();
            double[] resampledTestData = null;
            if (bTimesEqual) {
                // If time arrays for both sims are equal, no need to resample/interpolate, just obtain the datablock from dataManager
                SimDataBlock testSimDataBlock = testDataManager.getSimDataBlock(testVar.getName(), testTimeArray[j]);
                resampledTestData = testSimDataBlock.getData();
            } else {
                // Time resampling (interpolation) needed.
                while ((t < testTimeArray.length - 2) && (refTimeArray[j] >= testTimeArray[t + 1])) {
                    t++;
                }
                SimDataBlock testSimDataBlock_1 = testDataManager.getSimDataBlock(testVar.getName(), testTimeArray[t]);
                double[] testData_1 = testSimDataBlock_1.getData();
                SimDataBlock testSimDataBlock_2 = testDataManager.getSimDataBlock(testVar.getName(), testTimeArray[t + 1]);
                double[] testData_2 = testSimDataBlock_2.getData();
                resampledTestData = new double[testData_1.length];
                // 
                for (int m = 0; m < testData_1.length; m++) {
                    resampledTestData[m] = testData_1[m] + (testData_2[m] - testData_1[m]) * (refTimeArray[j] - testTimeArray[t]) / (testTimeArray[t + 1] - testTimeArray[t]);
                }
            }
            // Spatial resampling (interpolation) ...
            double[] spaceResampledData = new double[refData.length];
            if (!testMathDesc.getGeometry().getExtent().compareEqual(refMathDesc.getGeometry().getExtent()) || !testMathDesc.getGeometry().getOrigin().compareEqual(refMathDesc.getGeometry().getOrigin())) {
                throw new RuntimeException("Different origins and/or extents for the 2 geometries. Cannot compare the 2 simulations");
            }
            if (testMesh.getSizeX() != refMesh.getSizeX() || testMesh.getSizeY() != refMesh.getSizeY() || testMesh.getSizeZ() != refMesh.getSizeZ()) {
                if (testVar instanceof VolVariable) {
                    if (testMathDesc.getGeometry().getDimension() == 1 && refMathDesc.getGeometry().getDimension() == 1) {
                        spaceResampledData = resample1DSpatial(resampledTestData, testMesh, refMesh);
                    } else if (testMathDesc.getGeometry().getDimension() == 2 && refMathDesc.getGeometry().getDimension() == 2) {
                        spaceResampledData = resample2DSpatial(resampledTestData, testMesh, refMesh);
                    } else if (testMathDesc.getGeometry().getDimension() == 3 && refMathDesc.getGeometry().getDimension() == 3) {
                        spaceResampledData = resample3DSpatial(resampledTestData, testMesh, refMesh);
                    } else {
                        throw new RuntimeException("Comparison of 2 simulations with different geometry dimensions are not handled at this time!");
                    }
                } else {
                    throw new RuntimeException("spatial resampling for variable type: " + testVar.getClass().getName() + " not supported");
                }
            } else {
                // no space resampling required
                if (testVar instanceof MemVariable) {
                    // 
                    if (membraneIndexMapping == null) {
                        membraneIndexMapping = testMesh.getMembraneIndexMapping(refMesh);
                    }
                    spaceResampledData = new double[resampledTestData.length];
                    for (int k = 0; k < resampledTestData.length; k++) {
                        spaceResampledData[k] = resampledTestData[membraneIndexMapping[k]];
                    }
                } else {
                    // 
                    // no reordering needed for other variable types.
                    // 
                    spaceResampledData = resampledTestData;
                }
            }
            // for each point in data block ...
            testDataInfoProvider.getPDEDataContext().setVariableName(testVar.getName());
            for (int k = 0; k < refData.length; k++) {
                if (!testDataInfoProvider.isDefined(k)) {
                    continue;
                }
                // Determine maxRef, minRef, maxAbsErr for variable
                // SubDomain testSubDomain = null;
                String sn = null;
                VariableType refVarType = refDataID.getVariableType();
                if (refVarType.equals(VariableType.VOLUME)) {
                    // testSubDomain = refVolSubDomainLookup[k];
                    sn = refVolSubDomainLookup[k].getName();
                } else if (refVarType.equals(VariableType.MEMBRANE)) {
                    // testSubDomain = refMemSubDomainLookup[k];
                    sn = refMemSubDomainLookup[k].getName();
                } else if (refVarType.equals(VariableType.MEMBRANE_REGION)) {
                    sn = "MRV_" + i;
                } else if (refVarType.equals(VariableType.VOLUME_REGION)) {
                    sn = "VRV_" + i;
                } else {
                    throw new RuntimeException("Var " + refVar.getName() + " not supported yet!");
                }
                // hashKey = refVar.getName()+":"+testSubDomain.getName();
                hashKey = refVar.getName() + ":" + sn;
                tempVar = tempVarHash.get(hashKey);
                if (tempVar == null) {
                    tempVar = new DataErrorSummary(null);
                    tempVarHash.put(hashKey, tempVar);
                }
                tempVar.addDataValues(refData[k], spaceResampledData[k], refTimeArray[j], k, absErrorThreshold, relErrorThreshold);
            }
        // end for (k)
        }
    // end for (j)
    }
    // end for (i)
    Enumeration<String> enumKeys = tempVarHash.keys();
    while (enumKeys.hasMoreElements()) {
        String key = enumKeys.nextElement();
        DataErrorSummary tempVarSummary = tempVarHash.get(key);
        simComparisonSummary.addVariableComparisonSummary(new VariableComparisonSummary(key, tempVarSummary.getMinRef(), tempVarSummary.getMaxRef(), tempVarSummary.getMaxAbsoluteError(), tempVarSummary.getMaxRelativeError(), tempVarSummary.getL2Norm(), tempVarSummary.getTimeAtMaxAbsoluteError(), tempVarSummary.getIndexAtMaxAbsoluteError(), tempVarSummary.getTimeAtMaxRelativeError(), tempVarSummary.getIndexAtMaxRelativeError()));
    }
    return simComparisonSummary;
}

Example 33 with CompartmentSubDomain

use of cbit.vcell.math.CompartmentSubDomain in project vcell by virtualcell.

the class FiniteVolumeFileWriter method writeModelDescription.

/**
 *# Model description: FEATURE name handle priority boundary_conditions
 *MODEL_BEGIN
 *FEATURE Nucleus 0 value value value value
 *FEATURE Cytosol 1 flux flux value value
 *FEATURE ExtraCellular 2 flux flux value value
 *MEMBRANE Nucleus_Cytosol_membrane Nucleus Cytosol value value value value
 *MEMBRANE Cytosol_ExtraCellular_membrane Cytosol ExtraCellular flux flux value value
 *MODEL_END
 * @throws MathException
 */
private void writeModelDescription() throws MathException {
    Simulation simulation = simTask.getSimulation();
    printWriter.println("# Model description: FEATURE name handle boundary_conditions");
    printWriter.println(FVInputFileKeyword.MODEL_BEGIN);
    MathDescription mathDesc = simulation.getMathDescription();
    Enumeration<SubDomain> enum1 = mathDesc.getSubDomains();
    while (enum1.hasMoreElements()) {
        SubDomain sd = enum1.nextElement();
        if (sd instanceof CompartmentSubDomain) {
            CompartmentSubDomain csd = (CompartmentSubDomain) sd;
            printWriter.print(FVInputFileKeyword.FEATURE + " " + csd.getName() + " " + mathDesc.getHandle(csd) + " ");
            writeFeature_boundaryConditions(csd);
        } else if (sd instanceof MembraneSubDomain) {
            MembraneSubDomain msd = (MembraneSubDomain) sd;
            printWriter.print(FVInputFileKeyword.MEMBRANE + " " + msd.getName() + " " + msd.getInsideCompartment().getName() + " " + msd.getOutsideCompartment().getName() + " ");
            writeMembrane_boundaryConditions(msd);
        }
    }
    printWriter.println(FVInputFileKeyword.MODEL_END);
    printWriter.println();
}

Example 34 with CompartmentSubDomain

use of cbit.vcell.math.CompartmentSubDomain in project vcell by virtualcell.

the class FiniteVolumeFileWriter method writeMembrane_jumpConditions.

/**
 *JUMP_CONDITION_BEGIN r
 *INFLUX 0.0;
 *OUTFLUX 0.0;
 *JUMP_CONDITION_END
 * @throws ExpressionException
 * @throws MathException
 */
private void writeMembrane_jumpConditions(MembraneSubDomain msd) throws ExpressionException, MathException {
    Enumeration<JumpCondition> enum1 = msd.getJumpConditions();
    // equations for boundaryValues for inner compartment subdomain
    CompartmentSubDomain innerCompSubDomain = msd.getInsideCompartment();
    Enumeration<Equation> innercompSubDomEqnsEnum = innerCompSubDomain.getEquations();
    while (innercompSubDomEqnsEnum.hasMoreElements()) {
        Equation eqn = innercompSubDomEqnsEnum.nextElement();
        if (eqn instanceof PdeEquation) {
            PdeEquation pdeEqn = (PdeEquation) eqn;
            BoundaryConditionValue boundaryValue = pdeEqn.getBoundaryConditionValue(msd.getName());
            if (boundaryValue != null) {
                // check if the type of BoundaryConditionSpec for this membraneSubDomain (msd) in the (inner) compartmentSubDomain is Flux; if not, it cannot be handled.
                BoundaryConditionSpec bcs = innerCompSubDomain.getBoundaryConditionSpec(msd.getName());
                if (bcs == null) {
                    throw new MathException("No Boundary type specified for '" + msd.getName() + "' in '" + innerCompSubDomain.getName() + "'.");
                }
                if (bcs != null && !bcs.getBoundaryConditionType().compareEqual(BoundaryConditionType.getNEUMANN()) && !bChomboSolver) {
                    throw new MathException("Boundary type '" + bcs.getBoundaryConditionType().boundaryTypeStringValue() + "' for compartmentSubDomain '" + innerCompSubDomain.getName() + "' not handled by the chosen solver. Expecting boundary condition of type 'Flux'.");
                }
                if (pdeEqn.getVariable().getDomain() == null || pdeEqn.getVariable().getDomain().getName().equals(msd.getInsideCompartment().getName())) {
                    printWriter.println("JUMP_CONDITION_BEGIN " + pdeEqn.getVariable().getName());
                    Expression flux = subsituteExpression(boundaryValue.getBoundaryConditionExpression(), VariableDomain.VARIABLEDOMAIN_MEMBRANE);
                    String infix = replaceVolumeVariable(getSimulationTask(), msd, flux);
                    printWriter.println(bcs.getBoundaryConditionType().boundaryTypeStringValue().toUpperCase() + " " + msd.getInsideCompartment().getName() + " " + infix + ";");
                    printWriter.println("JUMP_CONDITION_END");
                    printWriter.println();
                }
            }
        }
    }
    // equations for boundaryValues for outer compartment subdomain
    CompartmentSubDomain outerCompSubDomain = msd.getOutsideCompartment();
    Enumeration<Equation> outerCompSubDomEqnsEnum = outerCompSubDomain.getEquations();
    while (outerCompSubDomEqnsEnum.hasMoreElements()) {
        Equation eqn = outerCompSubDomEqnsEnum.nextElement();
        if (eqn instanceof PdeEquation) {
            PdeEquation pdeEqn = (PdeEquation) eqn;
            BoundaryConditionValue boundaryValue = pdeEqn.getBoundaryConditionValue(msd.getName());
            if (boundaryValue != null) {
                // check if the type of BoundaryConditionSpec for this membraneSubDomain (msd) in the (inner) compartmentSubDomain is Flux; if not, it cannot be handled.
                BoundaryConditionSpec bcs = outerCompSubDomain.getBoundaryConditionSpec(msd.getName());
                if (bcs != null && !bcs.getBoundaryConditionType().compareEqual(BoundaryConditionType.getNEUMANN()) && !bChomboSolver) {
                    throw new MathException("Boundary type '" + bcs.getBoundaryConditionType().boundaryTypeStringValue() + "' for compartmentSubDomain '" + outerCompSubDomain.getName() + "' not handled by the chosen solver. Expecting boundary condition of type 'Flux'.");
                }
                if (pdeEqn.getVariable().getDomain() == null || pdeEqn.getVariable().getDomain().getName().equals(msd.getOutsideCompartment().getName())) {
                    printWriter.println("JUMP_CONDITION_BEGIN " + pdeEqn.getVariable().getName());
                    Expression flux = subsituteExpression(boundaryValue.getBoundaryConditionExpression(), VariableDomain.VARIABLEDOMAIN_MEMBRANE);
                    String infix = replaceVolumeVariable(getSimulationTask(), msd, flux);
                    printWriter.println(bcs.getBoundaryConditionType().boundaryTypeStringValue().toUpperCase() + " " + msd.getOutsideCompartment().getName() + " " + infix + ";");
                    printWriter.println("JUMP_CONDITION_END");
                    printWriter.println();
                }
            }
        }
    }
    while (enum1.hasMoreElements()) {
        JumpCondition jc = enum1.nextElement();
        printWriter.println("JUMP_CONDITION_BEGIN " + jc.getVariable().getName());
        // influx
        if (jc.getVariable().getDomain() == null || jc.getVariable().getDomain().getName().equals(msd.getInsideCompartment().getName())) {
            Expression flux = subsituteExpression(jc.getInFluxExpression(), VariableDomain.VARIABLEDOMAIN_MEMBRANE);
            String infix = replaceVolumeVariable(getSimulationTask(), msd, flux);
            printWriter.println(BoundaryConditionType.NEUMANN_STRING.toUpperCase() + " " + msd.getInsideCompartment().getName() + " " + infix + ";");
        }
        if (jc.getVariable().getDomain() == null || jc.getVariable().getDomain().getName().equals(msd.getOutsideCompartment().getName())) {
            // outflux
            Expression flux = subsituteExpression(jc.getOutFluxExpression(), VariableDomain.VARIABLEDOMAIN_MEMBRANE);
            String infix = replaceVolumeVariable(simTask, msd, flux);
            printWriter.println(BoundaryConditionType.NEUMANN_STRING.toUpperCase() + " " + msd.getOutsideCompartment().getName() + " " + infix + ";");
        }
        printWriter.println("JUMP_CONDITION_END");
        printWriter.println();
    }
}

Example 35 with CompartmentSubDomain

use of cbit.vcell.math.CompartmentSubDomain in project vcell by virtualcell.

the class FiniteVolumeFileWriter method writeFastSystem.

/*private void writeDataProcessor() throws DataAccessException, IOException, MathException, DivideByZeroException, ExpressionException {
	Simulation simulation = simTask.getSimulation();
	DataProcessingInstructions dpi = simulation.getDataProcessingInstructions();
	if (dpi == null) {
		printWriter.println("DATA_PROCESSOR_BEGIN " + DataProcessingInstructions.ROI_TIME_SERIES);
		printWriter.println("DATA_PROCESSOR_END");
		printWriter.println();
		return;
	}
	
	FieldDataIdentifierSpec fdis = dpi.getSampleImageFieldData(simulation.getVersion().getOwner());	
	if (fdis == null) {
		throw new DataAccessException("Can't find sample image in data processing instructions");
	}
	String secondarySimDataDir = PropertyLoader.getProperty(PropertyLoader.secondarySimDataDirProperty, null);	
	DataSetControllerImpl dsci = new DataSetControllerImpl(new NullSessionLog(),null,userDirectory.getParentFile(),secondarySimDataDir == null ? null : new File(secondarySimDataDir));
	CartesianMesh origMesh = dsci.getMesh(fdis.getExternalDataIdentifier());
	SimDataBlock simDataBlock = dsci.getSimDataBlock(null,fdis.getExternalDataIdentifier(), fdis.getFieldFuncArgs().getVariableName(), fdis.getFieldFuncArgs().getTime().evaluateConstant());
	VariableType varType = fdis.getFieldFuncArgs().getVariableType();
	VariableType dataVarType = simDataBlock.getVariableType();
	if (!varType.equals(VariableType.UNKNOWN) && !varType.equals(dataVarType)) {
		throw new IllegalArgumentException("field function variable type (" + varType.getTypeName() + ") doesn't match real variable type (" + dataVarType.getTypeName() + ")");
	}
	double[] origData = simDataBlock.getData();	

	String filename = SimulationJob.createSimulationJobID(Simulation.createSimulationID(simulation.getKey()), simulationJob.getJobIndex()) + FieldDataIdentifierSpec.getDefaultFieldDataFileNameForSimulation(fdis.getFieldFuncArgs());
	
	File fdatFile = new File(userDirectory, filename);
	
	
	DataSet.writeNew(fdatFile,
			new String[] {fdis.getFieldFuncArgs().getVariableName()},
			new VariableType[]{simDataBlock.getVariableType()},
			new ISize(origMesh.getSizeX(),origMesh.getSizeY(),origMesh.getSizeZ()),
			new double[][]{origData});
	printWriter.println("DATA_PROCESSOR_BEGIN " + dpi.getScriptName());
	printWriter.println(dpi.getScriptInput());
	printWriter.println("SampleImageFile " + fdis.getFieldFuncArgs().getVariableName() + " " + fdis.getFieldFuncArgs().getTime().infix() + " " + fdatFile);
	printWriter.println("DATA_PROCESSOR_END");
	printWriter.println();

	
}
*/
/**
 *# fast system dimension num_dependents
 *FAST_SYSTEM_BEGIN 2 2
 *INDEPENDENT_VARIALBES rf r
 *DEPENDENT_VARIALBES rB rfB
 *
 *PSEUDO_CONSTANT_BEGIN
 *__C0 (rfB + rf);
 *__C1 (r + rB);
 *PSEUDO_CONSTANT_END
 *
 *FAST_RATE_BEGIN
 * - ((0.02 * ( - ( - r + __C1) - ( - rf + __C0) + (20.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))) + _VCell_FieldData_0) * rf) - (0.1 * ( - rf + __C0)));
 *((0.02 * r * ( - ( - r + __C1) - ( - rf + __C0) + (20.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))) + _VCell_FieldData_0)) - (0.1 * ( - r + __C1)));
 *FAST_RATE_END
 *
 *FAST_DEPENDENCY_BEGIN
 *rB ( - r + __C1);
 *rfB ( - rf + __C0);
 *FAST_DEPENDENCY_END
 *
 *JACOBIAN_BEGIN
 * - (0.1 + (0.02 * (1.0 + (0.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0)))) * rf) + (0.02 * ( - ( - r + __C1) - ( - rf + __C0) + (20.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))) + _VCell_FieldData_0)));
 * - (0.02 * (1.0 + (0.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0)))) * rf);
 *(0.02 * r * (1.0 + (0.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0)))));
 *(0.1 + (0.02 * ( - ( - r + __C1) - ( - rf + __C0) + (20.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))) + _VCell_FieldData_0)) + (0.02 * r * (1.0 + (0.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))))));
 *JACOBIAN_END
 *
 *FAST_SYSTEM_END
 * @throws ExpressionException
 * @throws MathException
 */
private void writeFastSystem(SubDomain subDomain) throws MathException, ExpressionException {
    VariableDomain variableDomain = (subDomain instanceof CompartmentSubDomain) ? VariableDomain.VARIABLEDOMAIN_VOLUME : VariableDomain.VARIABLEDOMAIN_MEMBRANE;
    FastSystem fastSystem = subDomain.getFastSystem();
    if (fastSystem == null) {
        return;
    }
    FastSystemAnalyzer fs_analyzer = new FastSystemAnalyzer(fastSystem, simTask.getSimulationJob().getSimulationSymbolTable());
    int numIndep = fs_analyzer.getNumIndependentVariables();
    int numDep = fs_analyzer.getNumDependentVariables();
    int numPseudo = fs_analyzer.getNumPseudoConstants();
    printWriter.println("# fast system dimension num_dependents");
    printWriter.println("FAST_SYSTEM_BEGIN " + numIndep + " " + numDep);
    if (numIndep != 0) {
        printWriter.print("INDEPENDENT_VARIALBES ");
        Enumeration<Variable> enum1 = fs_analyzer.getIndependentVariables();
        while (enum1.hasMoreElements()) {
            Variable var = enum1.nextElement();
            printWriter.print(var.getName() + " ");
        }
        printWriter.println();
    }
    if (numDep != 0) {
        printWriter.print("DEPENDENT_VARIALBES ");
        Enumeration<Variable> enum1 = fs_analyzer.getDependentVariables();
        while (enum1.hasMoreElements()) {
            Variable var = enum1.nextElement();
            printWriter.print(var.getName() + " ");
        }
        printWriter.println();
    }
    printWriter.println();
    if (numPseudo != 0) {
        printWriter.println("PSEUDO_CONSTANT_BEGIN");
        Enumeration<PseudoConstant> enum1 = fs_analyzer.getPseudoConstants();
        while (enum1.hasMoreElements()) {
            PseudoConstant pc = enum1.nextElement();
            printWriter.println(pc.getName() + " " + subsituteExpression(pc.getPseudoExpression(), fs_analyzer, variableDomain).infix() + ";");
        }
        printWriter.println("PSEUDO_CONSTANT_END");
        printWriter.println();
    }
    if (numIndep != 0) {
        printWriter.println("FAST_RATE_BEGIN");
        Enumeration<Expression> enum1 = fs_analyzer.getFastRateExpressions();
        while (enum1.hasMoreElements()) {
            Expression exp = enum1.nextElement();
            printWriter.println(subsituteExpression(exp, fs_analyzer, variableDomain).infix() + ";");
        }
        printWriter.println("FAST_RATE_END");
        printWriter.println();
    }
    if (numDep != 0) {
        printWriter.println("FAST_DEPENDENCY_BEGIN");
        Enumeration<Expression> enum_exp = fs_analyzer.getDependencyExps();
        Enumeration<Variable> enum_var = fs_analyzer.getDependentVariables();
        while (enum_exp.hasMoreElements()) {
            Expression exp = enum_exp.nextElement();
            Variable depVar = enum_var.nextElement();
            printWriter.println(depVar.getName() + " " + subsituteExpression(exp, fs_analyzer, variableDomain).infix() + ";");
        }
        printWriter.println("FAST_DEPENDENCY_END");
        printWriter.println();
    }
    if (numIndep != 0) {
        printWriter.println("JACOBIAN_BEGIN");
        Enumeration<Expression> enum_fre = fs_analyzer.getFastRateExpressions();
        while (enum_fre.hasMoreElements()) {
            Expression fre = enum_fre.nextElement();
            Enumeration<Variable> enum_var = fs_analyzer.getIndependentVariables();
            while (enum_var.hasMoreElements()) {
                Variable var = enum_var.nextElement();
                Expression exp = subsituteExpression(fre, fs_analyzer, variableDomain).flatten();
                Expression differential = exp.differentiate(var.getName());
                printWriter.println(subsituteExpression(differential, fs_analyzer, variableDomain).infix() + ";");
            }
        }
        printWriter.println("JACOBIAN_END");
        printWriter.println();
    }
    printWriter.println("FAST_SYSTEM_END");
    printWriter.println();
}