Search in sources :

Example 1 with VCCDifference

use of org.vcell.solver.comsol.model.VCCGeomFeature.VCCDifference in project vcell by virtualcell.

the class ComsolModelBuilder method csgVisitor.

private static VCCGeomFeature csgVisitor(CSGNode node, ArrayList<VCCGeomFeature> geomFeatureList, String namePrefix) {
    final VCCGeomFeature newFeature;
    if (node instanceof CSGPrimitive) {
        CSGPrimitive csg = (CSGPrimitive) node;
        switch(csg.getType()) {
            case CONE:
                {
                    newFeature = new VCCCone(namePrefix + csg.getName());
                    break;
                }
            case CUBE:
                {
                    String[] size = new String[] { "2", "2", "2" };
                    String[] pos = new String[] { "-1", "-1", "-1" };
                    VCCBlock vccblock = new VCCBlock(namePrefix + csg.getName(), size, pos);
                    newFeature = vccblock;
                    break;
                }
            case CYLINDER:
                {
                    newFeature = new VCCCylinder(namePrefix + csg.getName());
                    break;
                }
            case SPHERE:
                {
                    String r = "1";
                    String[] pos = new String[] { "0", "0", "0" };
                    newFeature = new VCCSphere(namePrefix + csg.getName(), pos, r);
                    break;
                }
            default:
                {
                    throw new RuntimeException("csg primative type '" + csg.getType().name() + "' not yet supported in COMSOL model builder");
                }
        }
    } else if (node instanceof CSGSetOperator) {
        CSGSetOperator setOp = (CSGSetOperator) node;
        ArrayList<VCCGeomFeature> childSubtrees = new ArrayList<VCCGeomFeature>();
        for (CSGNode child : setOp.getChildren()) {
            VCCGeomFeature vccChild = csgVisitor(child, geomFeatureList, namePrefix);
            childSubtrees.add(vccChild);
        }
        switch(setOp.getOpType()) {
            case DIFFERENCE:
                {
                    if (childSubtrees.size() != 2) {
                        throw new RuntimeException("expecting exactly two children for CSG difference operator");
                    }
                    VCCDifference diff = new VCCDifference(namePrefix + setOp.getName(), Keep.off);
                    diff.input.add(childSubtrees.get(0));
                    diff.input2.add(childSubtrees.get(1));
                    newFeature = diff;
                    break;
                }
            case INTERSECTION:
                {
                    if (childSubtrees.size() < 2) {
                        throw new RuntimeException("expecting two or more children for CSG intersection operator");
                    }
                    VCCIntersection intersection = new VCCIntersection(namePrefix + setOp.getName(), Keep.off);
                    intersection.input.add(childSubtrees.get(0));
                    newFeature = intersection;
                    break;
                }
            case UNION:
                {
                    if (childSubtrees.size() < 2) {
                        throw new RuntimeException("expecting two or more children for CSG union operator");
                    }
                    VCCUnion union = new VCCUnion(namePrefix + setOp.getName(), Keep.off);
                    union.input.add(childSubtrees.get(0));
                    newFeature = union;
                    break;
                }
            default:
                {
                    throw new RuntimeException("csg set operator '" + setOp.getOpType().name() + "' not yet supported in COMSOL model builder");
                }
        }
    } else if (node instanceof CSGTransformation) {
        CSGTransformation transformation = (CSGTransformation) node;
        VCCGeomFeature vccChild = csgVisitor(transformation.getChild(), geomFeatureList, namePrefix);
        if (transformation instanceof CSGHomogeneousTransformation) {
            throw new RuntimeException("unsupported CSG transformation type Homogeneous transformation");
        } else if (transformation instanceof CSGRotation) {
            CSGRotation rotation = (CSGRotation) transformation;
            String[] axis = new String[] { Double.toString(rotation.getAxis().getX()), Double.toString(rotation.getAxis().getY()), Double.toString(rotation.getAxis().getZ()) };
            String[] pos = new String[] { "0.0", "0.0", "0.0" };
            String rot = Double.toString(rotation.getRotationRadians());
            VCCRotate vccrotate = new VCCRotate(namePrefix + rotation.getName(), axis, pos, rot, Keep.off);
            vccrotate.input.add(vccChild);
            newFeature = vccrotate;
        } else if (transformation instanceof CSGScale) {
            CSGScale scale = (CSGScale) transformation;
            String[] factor = new String[] { Double.toString(scale.getScale().getX()), Double.toString(scale.getScale().getY()), Double.toString(scale.getScale().getZ()) };
            String[] pos = new String[] { "0.0", "0.0", "0.0" };
            VCCScale vccscale = new VCCScale(namePrefix + scale.getName(), pos, factor, Keep.off);
            vccscale.input.add(vccChild);
            newFeature = vccscale;
        } else if (transformation instanceof CSGTranslation) {
            CSGTranslation translation = (CSGTranslation) transformation;
            String[] displ = new String[] { Double.toString(translation.getTranslation().getX()), Double.toString(translation.getTranslation().getY()), Double.toString(translation.getTranslation().getZ()) };
            VCCMove vccmove = new VCCMove(namePrefix + translation.getName(), displ, Keep.off);
            vccmove.input.add(vccChild);
            newFeature = vccmove;
        } else {
            throw new RuntimeException("unsupported CSG transformation type '" + transformation.getClass().getSimpleName() + "'");
        }
    } else {
        throw new RuntimeException("unsupported CSGNode type '" + node.getClass().getSimpleName() + "'");
    }
    geomFeatureList.add(newFeature);
    return newFeature;
}
Also used : VCCDifference(org.vcell.solver.comsol.model.VCCGeomFeature.VCCDifference) VCCCone(org.vcell.solver.comsol.model.VCCGeomFeature.VCCCone) VCCUnion(org.vcell.solver.comsol.model.VCCGeomFeature.VCCUnion) CSGTranslation(cbit.vcell.geometry.CSGTranslation) CSGTransformation(cbit.vcell.geometry.CSGTransformation) VCCMove(org.vcell.solver.comsol.model.VCCGeomFeature.VCCMove) CSGPrimitive(cbit.vcell.geometry.CSGPrimitive) VCCBlock(org.vcell.solver.comsol.model.VCCGeomFeature.VCCBlock) ArrayList(java.util.ArrayList) CSGNode(cbit.vcell.geometry.CSGNode) VCCSphere(org.vcell.solver.comsol.model.VCCGeomFeature.VCCSphere) VCCIntersection(org.vcell.solver.comsol.model.VCCGeomFeature.VCCIntersection) CSGRotation(cbit.vcell.geometry.CSGRotation) CSGScale(cbit.vcell.geometry.CSGScale) VCCScale(org.vcell.solver.comsol.model.VCCGeomFeature.VCCScale) VCCCylinder(org.vcell.solver.comsol.model.VCCGeomFeature.VCCCylinder) CSGHomogeneousTransformation(cbit.vcell.geometry.CSGHomogeneousTransformation) VCCRotate(org.vcell.solver.comsol.model.VCCGeomFeature.VCCRotate) VCCGeomFeature(org.vcell.solver.comsol.model.VCCGeomFeature) CSGSetOperator(cbit.vcell.geometry.CSGSetOperator)

Example 2 with VCCDifference

use of org.vcell.solver.comsol.model.VCCGeomFeature.VCCDifference in project vcell by virtualcell.

the class ComsolServiceScripting method run.

// private String tojs(String[] a){
// StringBuffer buffer = new StringBuffer();
// buffer.append("[");
// for (int i=0;i<a.length;i++){
// buffer.append("'"+a[i]+"'");
// if (i<a.length-1){
// buffer.append(",");
// }
// }
// buffer.append("]");
// return buffer.toString();
// }
private void run(ScriptEngine engine, VCCModel vccModel, File reportFile, File javaFile, File mphFile, File jsFile) throws ScriptException, IOException {
    StringBuffer buffer = new StringBuffer();
    // ModelUtil.initStandalone(true);
    buffer.append("com.comsol.model.util.ModelUtil.initStandalone(true);" + "\n");
    // Model model =                       ModelUtil.create(vccModel.name);
    buffer.append("model = com.comsol.model.util.ModelUtil.create('" + vccModel.name + "');" + "\n");
    // //               model.modelPath(vccModel.modelpath);
    // buffer.append(	"model.modelPath('"+vccModel.modelpath+"');"	+ "\n");
    // model.comments(vccModel.comments);
    buffer.append("model.comments('" + vccModel.comments.replaceAll("\n", " ") + "');" + "\n");
    for (VCCModelNode modelNode : vccModel.modelnodes) {
        // model.modelNode().create(modelNode.name);
        buffer.append("model.modelNode().create('" + modelNode.name + "');" + "\n");
    }
    boolean bBeforeSelections = true;
    for (VCCGeomSequence geomSequence : vccModel.geometrysequences) {
        // model.geom().create(geomSequence.name, geomSequence.dim);
        buffer.append("model.geom().create('" + geomSequence.name + "'," + geomSequence.dim + ");" + "\n");
        for (VCCGeomFeature geomFeature : geomSequence.geomfeatures) {
            switch(geomFeature.type) {
                case Circle:
                    {
                        buffer.append("print('geomFeature Circle named " + geomFeature.name + "');" + "\n");
                        VCCCircle circle = (VCCCircle) geomFeature;
                        // model.geom(geomSequence.name)      .create(circle.name, "Circle");
                        buffer.append("model.geom('" + geomSequence.name + "').create('" + circle.name + "','Circle');" + "\n");
                        // model.geom(geomSequence.name)      .feature(circle.name)      .set("selresult", "on");
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + circle.name + "').set('selresult','on');" + "\n");
                        // model.geom(geomSequence.name)      .feature(circle.name)      .set("pos", circle.pos);
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + circle.name + "').set('pos', " + tojs(circle.pos) + ");" + "\n");
                        // model.geom(geomSequence.name)		.feature(circle.name)	   .set("r", circle.r);
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + circle.name + "').set('r', '" + circle.r + "');" + "\n");
                        break;
                    }
                case Sphere:
                    {
                        buffer.append("print('geomFeature Sphere named " + geomFeature.name + "');" + "\n");
                        VCCSphere sphere = (VCCSphere) geomFeature;
                        // model.geom(geomSequence.name)      .create(sphere.name, "Sphere");
                        buffer.append("model.geom('" + geomSequence.name + "').create('" + sphere.name + "','Sphere');" + "\n");
                        // model.geom(geomSequence.name)      .feature(sphere.name)      .set("pos", sphere.pos);
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + sphere.name + "').set('pos', " + tojs(sphere.pos) + ");" + "\n");
                        // model.geom(geomSequence.name)      .feature(sphere.name)      .set("r", sphere.r);
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + sphere.name + "').set('r', " + sphere.r + ");" + "\n");
                        break;
                    }
                case Square:
                    {
                        VCCSquare square = (VCCSquare) geomFeature;
                        // model.geom(geomSequence.name)      .create(square.name, "Square");
                        buffer.append("model.geom('" + geomSequence.name + "').create('" + square.name + "','Square');" + "\n");
                        // buffer.append(	"model.geom('"+geomSequence.name+"').feature('"+square.name+"').set('selresult','on');"	+ "\n");
                        break;
                    }
                case Block:
                    {
                        buffer.append("print('geomFeature Block named " + geomFeature.name + "');" + "\n");
                        VCCBlock block = (VCCBlock) geomFeature;
                        // model.geom(geomSequence.name)      .create(block.name, "Block");
                        buffer.append("model.geom('" + geomSequence.name + "').create('" + block.name + "','Block');" + "\n");
                        // model.geom(geomSequence.name)      .feature(block.name)      .set("pos", block.pos);
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + block.name + "').set('pos', " + tojs(block.pos) + ");" + "\n");
                        // model.geom(geomSequence.name)      .feature(block.name)      .set("size", bock.size);
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + block.name + "').set('size', " + tojs(block.size) + ");" + "\n");
                        break;
                    }
                case Difference:
                    {
                        buffer.append("print('geomFeature Difference named " + geomFeature.name + "');" + "\n");
                        VCCDifference diff = (VCCDifference) geomFeature;
                        // model.geom(geomSequence.name)      .create(diff.name, "Difference");
                        buffer.append("model.geom('" + geomSequence.name + "').create('" + diff.name + "','Difference');" + "\n");
                        // model.geom(geomSequence.name)      .feature(diff.name)      .set("selresult", "on");
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + diff.name + "').set('selresult','on');" + "\n");
                        ArrayList<String> inputNames = new ArrayList<String>();
                        for (VCCGeomFeature feature : diff.input) {
                            inputNames.add(feature.name);
                        }
                        String[] inputSet = inputNames.toArray(new String[0]);
                        ArrayList<String> input2Names = new ArrayList<String>();
                        for (VCCGeomFeature feature : diff.input2) {
                            input2Names.add(feature.name);
                        }
                        String[] input2Set = input2Names.toArray(new String[0]);
                        // model.geom(geomSequence.name)      .feature(diff.name)      .selection("input").set(inputSet);
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + diff.name + "').selection('input').set(" + tojs(inputSet) + ");" + "\n");
                        // model.geom(geomSequence.name)      .feature(diff.name)      .selection("input2").set(input2Set);
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + diff.name + "').selection('input2').set(" + tojs(input2Set) + ");" + "\n");
                        if (diff.keep == Keep.on) {
                            // model.geom(geomSequence.name)      .feature(diff.name)      .set("keep", "on");
                            buffer.append("model.geom('" + geomSequence.name + "').feature('" + diff.name + "').set('keep', 'on');" + "\n");
                        } else if (diff.keep == Keep.off) {
                            // model.geom(geomSequence.name)      .feature(diff.name)      .set("keep", "off");
                            buffer.append("model.geom('" + geomSequence.name + "').feature('" + diff.name + "').set('keep', 'off');" + "\n");
                        }
                        break;
                    }
                case Scale:
                    {
                        buffer.append("print('geomFeature Scale named " + geomFeature.name + "');" + "\n");
                        VCCScale scale = (VCCScale) geomFeature;
                        // model.geom(geomSequence.name)      .create(diff.name, "Difference");
                        buffer.append("model.geom('" + geomSequence.name + "').create('" + scale.name + "','Scale');" + "\n");
                        ArrayList<String> inputNames = new ArrayList<String>();
                        for (VCCGeomFeature feature : scale.input) {
                            inputNames.add(feature.name);
                        }
                        String[] inputSet = inputNames.toArray(new String[0]);
                        // model.geom(geomSequence.name)      .feature(scale.name)      .selection("input").set(inputSet);
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + scale.name + "').selection('input').set(" + tojs(inputSet) + ");" + "\n");
                        // model.geom(geomSequence.name)      .feature(scale.name)      .set('factor', scale.factor);
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + scale.name + "').set('factor'," + tojs(scale.factor) + ");" + "\n");
                        break;
                    }
                case Move:
                    {
                        buffer.append("print('geomFeature Move named " + geomFeature.name + "');" + "\n");
                        VCCMove move = (VCCMove) geomFeature;
                        // model.geom(geomSequence.name)      .create(diff.name, "Difference");
                        buffer.append("model.geom('" + geomSequence.name + "').create('" + move.name + "','Move');" + "\n");
                        ArrayList<String> inputNames = new ArrayList<String>();
                        for (VCCGeomFeature feature : move.input) {
                            inputNames.add(feature.name);
                        }
                        String[] inputSet = inputNames.toArray(new String[0]);
                        // model.geom(geomSequence.name)      .feature(move.name)      .selection("input").set(inputSet);
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + move.name + "').selection('input').set(" + tojs(inputSet) + ");" + "\n");
                        // model.geom(geomSequence.name)      .feature(move.name)      .set("displ", move.displ);
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + move.name + "').set('displ'," + tojs(move.displ) + ");" + "\n");
                        break;
                    }
                case IntersectionSelection:
                    {
                        buffer.append("print('geomFeature IntersectionSelection named " + geomFeature.name + "');" + "\n");
                        if (bBeforeSelections) {
                            bBeforeSelections = false;
                            // model.geom("geom1")                .run("fin");
                            buffer.append("model.geom('" + geomSequence.name + "').run('fin');" + "\n");
                        }
                        VCCIntersectionSelection intersectionSelection = (VCCIntersectionSelection) geomFeature;
                        // model.geom("geom1")                .create("intsel1", "IntersectionSelection");
                        buffer.append("model.geom('" + geomSequence.name + "').create('" + intersectionSelection.name + "', 'IntersectionSelection');" + "\n");
                        // model.geom("geom1")                .feature("intsel1")                       .set("entitydim", "2");
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + intersectionSelection.name + "').set('entitydim', '" + intersectionSelection.entitydim + "');" + "\n");
                        ArrayList<String> inputNames = new ArrayList<String>();
                        for (VCCGeomFeature feature : intersectionSelection.input) {
                            inputNames.add(feature.name);
                        }
                        String[] inputSet = inputNames.toArray(new String[0]);
                        // model.geom("geom1")                .feature("intsel1")                       .set("input", new String[]{"diffec", "celltranslation0"});
                        buffer.append("model.geom('" + geomSequence.name + "').feature('" + intersectionSelection.name + "').set('input'," + tojs(inputSet) + ");" + "\n");
                        break;
                    }
                default:
                    {
                        throw new RuntimeException("unexpected Geometry Feature type " + geomFeature.type.name());
                    }
            }
            if (!(geomFeature instanceof VCCIntersectionSelection)) {
                // model.geom(geomSequence.name)		 .feature(geomFeature.name)		 .set("selresult", "on");
                buffer.append("model.geom('" + geomSequence.name + "').feature('" + geomFeature.name + "').set('selresult','on');" + "\n");
            }
        }
        System.out.println("building geometry");
        // long t_startGeomRun_ms = System.currentTimeMillis();
        buffer.append("var t_startGeomRun_ms = Date.now();" + "\n");
        // model.geom(geomSequence.name)      .run();
        buffer.append("model.geom('" + geomSequence.name + "').run();" + "\n");
        // long t_endGeomRun_ms = System.currentTimeMillis();
        buffer.append("var t_endGeomRun_ms = Date.now();" + "\n");
        // System.out.println("geometry built in " + ((t_endGeomRun_ms - t_startGeomRun_ms) / 1000.0) + " seconds");
        buffer.append("print('geometry built in ' + ((t_endGeomRun_ms - t_startGeomRun_ms) / 1000.0) + ' seconds');" + "\n");
    }
    for (VCCMeshSequence meshSequence : vccModel.meshes) {
        // model.mesh().create(meshSequence.name, meshSequence.geom.name);
        buffer.append("model.mesh().create('" + meshSequence.name + "','" + meshSequence.geom.name + "');" + "\n");
    }
    for (VCCPhysics physics : vccModel.physics) {
        switch(physics.type) {
            // }
            case ConvectionDiffusionEquation:
                {
                    VCCConvectionDiffusionEquation pde = (VCCConvectionDiffusionEquation) physics;
                    buffer.append("print('(0) physics ConvectionDiffusionEquation for field name " + pde.fieldName + "');" + "\n");
                    // model.physics().create(pde.name, "ConvectionDiffusionEquation", pde.geom.name);
                    buffer.append("model.physics().create('" + pde.name + "','ConvectionDiffusionEquation','" + pde.geom.name + "');" + "\n");
                    buffer.append("print('(1) physics ConvectionDiffusionEquation for field name " + pde.fieldName + "');" + "\n");
                    String c = pde.diffTerm_c;
                    if (c != null) {
                        if (vccModel.dim == 2) {
                            // model.physics(pde.name)      .feature("cdeq1").setIndex("c", new String[] { c, "0", "0", c }, 0);
                            buffer.append("model.physics('" + pde.name + "').feature('cdeq1').setIndex('c', " + tojs(new String[] { c, "0", "0", c }) + ", 0);" + "\n");
                        } else if (vccModel.dim == 3) {
                            // model.physics(pde.name)      .feature("cdeq1").setIndex("c", new String[] { c, "0", "0", c }, 0);
                            buffer.append("model.physics('" + pde.name + "').feature('cdeq1').setIndex('c', " + tojs(new String[] { c, "0", "0", "0", c, "0", "0", "0", c }) + ", 0);" + "\n");
                        }
                    }
                    if (pde.advection_be != null) {
                        // model.physics(pde.name)      .feature("cdeq1").setIndex("be", pde.advection_be, 0);
                        buffer.append("model.physics('" + pde.name + "').feature('cdeq1').setIndex('be', " + tojs(pde.advection_be) + ", 0);" + "\n");
                    }
                    buffer.append("print('(2) physics ConvectionDiffusionEquation for field name " + pde.fieldName + "');" + "\n");
                    // model.physics(pde.name)      .feature("cdeq1").setIndex("f", pde.sourceTerm_f, 0);
                    buffer.append("model.physics('" + pde.name + "').feature('cdeq1').setIndex('f', '" + pde.sourceTerm_f + "', 0);" + "\n");
                    buffer.append("print('(3) physics ConvectionDiffusionEquation for field name " + pde.fieldName + "');" + "\n");
                    // model.physics(pde.name)      .field("dimensionless")  .field(pde.fieldName);
                    buffer.append("model.physics('" + pde.name + "').field('dimensionless').field('" + pde.fieldName + "');" + "\n");
                    buffer.append("print('(4) physics ConvectionDiffusionEquation for field name " + pde.fieldName + "');" + "\n");
                    // model.physics(pde.name)      .feature("init1").set(pde.fieldName, pde.initial);
                    buffer.append("model.physics('" + pde.name + "').feature('init1').set('" + pde.fieldName + "','" + pde.initial + "');" + "\n");
                    String selectionName = null;
                    if (pde.dim == vccModel.dim) {
                        selectionName = pde.geom.name + "_" + pde.geomFeature.name + "_dom";
                    } else if (pde.dim == vccModel.dim - 1) {
                        selectionName = pde.geom.name + "_" + pde.geomFeature.name + "_bnd";
                    }
                    buffer.append("print('(5) physics ConvectionDiffusionEquation for field name " + pde.fieldName + ", selection " + selectionName + "');" + "\n");
                    // model.physics(pde.name)      .selection().named(selectionName);
                    buffer.append("model.physics('" + pde.name + "').selection().named('" + selectionName + "');" + "\n");
                    buffer.append("print('(6) physics ConvectionDiffusionEquation for field name " + pde.fieldName + ", selection " + selectionName + "');" + "\n");
                    for (VCCPhysicsFeature feature : pde.features) {
                        if (feature instanceof VCCFluxBoundary) {
                            VCCFluxBoundary fluxBoundary = (VCCFluxBoundary) feature;
                            String boundarySelectionName = pde.geom.name + "_" + fluxBoundary.selection.name;
                            buffer.append("print('(0) physics ConvectionDiffusionEquation for field name " + pde.fieldName + ", FluxBoundary name " + fluxBoundary.name + ", selection name " + boundarySelectionName + "');" + "\n");
                            // model.physics(pde.name)      .create(fluxBoundary.name, "FluxBoundary", 2);
                            buffer.append("model.physics('" + pde.name + "').create('" + fluxBoundary.name + "','FluxBoundary'," + fluxBoundary.dim + ");" + "\n");
                            buffer.append("print('(1) physics ConvectionDiffusionEquation for field name " + pde.fieldName + ", FluxBoundary name " + fluxBoundary.name + ", selection name " + boundarySelectionName + "');" + "\n");
                            // model.physics(pde.name)      .feature(fluxBoundary.name)      .selection().named("geom1_intsel1");
                            buffer.append("model.physics('" + pde.name + "').feature('" + fluxBoundary.name + "').selection().named('" + boundarySelectionName + "');" + "\n");
                            buffer.append("print('(2) physics ConvectionDiffusionEquation for field name " + pde.fieldName + ", FluxBoundary name " + fluxBoundary.name + ", selection name " + boundarySelectionName + "');" + "\n");
                            // model.physics(pde.name)      .feature(fluxBoundary.name)      .set("g", "5-RanC_nuc");
                            buffer.append("model.physics('" + pde.name + "').feature('" + fluxBoundary.name + "').set('g', '" + fluxBoundary.flux_g + "');" + "\n");
                            buffer.append("print('(3) physics ConvectionDiffusionEquation for field name " + pde.fieldName + ", FluxBoundary name " + fluxBoundary.name + ", selection name " + boundarySelectionName + "');" + "\n");
                        }
                    }
                    break;
                }
            default:
                {
                    throw new RuntimeException("unsupported physics type " + physics.type.name() + " in COMSOL script builder");
                }
        }
    }
    System.out.println("building mesh");
    // long t_startMeshRun_ms = System.currentTimeMillis();
    buffer.append("var t_startMeshRun_ms = Date.now();" + "\n");
    // model.mesh("mesh1").run();
    buffer.append("model.mesh('mesh1').run();" + "\n");
    // long t_endMeshRun_ms = System.currentTimeMillis();
    buffer.append("var t_endMeshRun_ms = Date.now();" + "\n");
    // System.out.println("mesh built in " + ((t_endMeshRun_ms - t_startMeshRun_ms) / 1000.0) + " seconds");
    buffer.append("print('mesh built in ' + ((t_endMeshRun_ms - t_startMeshRun_ms) / 1000.0) + ' seconds');" + "\n");
    // model.study().create(vccModel.study.name);
    buffer.append("model.study().create('" + vccModel.study.name + "');" + "\n");
    for (VCCStudyFeature studyFeature : vccModel.study.features) {
        switch(studyFeature.type) {
            case Transient:
                {
                    VCCTransientStudyFeature transientTime = (VCCTransientStudyFeature) studyFeature;
                    // model.study(vccModel.study.name)      .create("time", "Transient");
                    buffer.append("model.study('" + vccModel.study.name + "').create('time', 'Transient');" + "\n");
                    String range = "range(" + transientTime.startingTime + "," + transientTime.timeStep + "," + transientTime.endTime + ")";
                    // model.study(vccModel.study.name)      .feature("time").set("tlist", range);
                    buffer.append("model.study('" + vccModel.study.name + "').feature('time').set('tlist', '" + range + "');" + "\n");
                    break;
                }
            default:
                {
                    throw new RuntimeException("unknown studyFeature type " + studyFeature.type.name());
                }
        }
    // model.study(study.name).create(studyFeature.name,
    // studyFeature.type.name());
    // for (VCCPhysics activePhysics : studyFeature.activePhysics){
    // model.study(study.name).feature(studyFeature.name).activate(activePhysics.name);
    // }
    }
    // model.sol().create("sol1");
    buffer.append("model.sol().create('sol1');" + "\n");
    // model.sol("sol1").study(vccModel.study.name);
    buffer.append("model.sol('sol1').study('" + vccModel.study.name + "');" + "\n");
    // model.sol("sol1").attach(vccModel.study.name);
    buffer.append("model.sol('sol1').attach('" + vccModel.study.name + "');" + "\n");
    // model.sol("sol1").create("st1", "StudyStep");
    buffer.append("model.sol('sol1').create('st1', 'StudyStep');" + "\n");
    // model.sol("sol1").create("v1", "Variables");
    buffer.append("model.sol('sol1').create('v1', 'Variables');" + "\n");
    // model.sol("sol1").create("t1", "Time");
    buffer.append("model.sol('sol1').create('t1', 'Time');" + "\n");
    // model.sol("sol1").feature("t1").create("fc1", "FullyCoupled");
    buffer.append("model.sol('sol1').feature('t1').create('fc1', 'FullyCoupled');" + "\n");
    // model.sol("sol1").feature("t1").feature().remove("fcDef");
    buffer.append("model.sol('sol1').feature('t1').feature().remove('fcDef');" + "\n");
    // model.result().export().create("data1", "Data");
    buffer.append("model.result().export().create('data1', 'Data');" + "\n");
    // model.sol("sol1").attach(vccModel.study.name);
    buffer.append("model.sol('sol1').attach('" + vccModel.study.name + "');" + "\n");
    // long t_beginsolver_ms = System.currentTimeMillis();
    buffer.append("var t_beginsolver_ms = Date.now();" + "\n");
    // System.out.println("running solver");
    buffer.append("print('running solver');" + "\n");
    // model.sol("sol1").runAll();
    buffer.append("model.sol('sol1').runAll();" + "\n");
    // long t_endsolver_ms = System.currentTimeMillis();
    buffer.append("var t_endsolver_ms = Date.now();" + "\n");
    // System.out.println("solver finished in " + ((t_endsolver_ms - t_beginsolver_ms) / 1000.0) + " seconds");
    buffer.append("print('solver finished in ' + ((t_endsolver_ms - t_beginsolver_ms) / 1000.0) + ' seconds');" + "\n");
    // model.result().export("data1").label("Data 2");
    buffer.append("model.result().export('data1').label('Data 2');" + "\n");
    // model.result().export("data1").set("struct", "sectionwise");
    buffer.append("model.result().export('data1').set('struct', 'sectionwise');" + "\n");
    // model.result().export("data1").set("filename", reportFile.getAbsolutePath());
    String reportFileName = StringEscapeUtils.escapeEcmaScript(reportFile.getAbsolutePath());
    buffer.append("model.result().export('data1').set('filename', '" + reportFileName + "');" + "\n");
    ArrayList<String> units = new ArrayList<String>();
    ArrayList<String> descriptions = new ArrayList<String>();
    ArrayList<String> expressions = new ArrayList<String>();
    for (VCCPhysics physics : vccModel.physics) {
        if (physics instanceof VCCConvectionDiffusionEquation) {
            VCCConvectionDiffusionEquation pde = (VCCConvectionDiffusionEquation) physics;
            units.add("1");
            descriptions.add(pde.fieldName);
            expressions.add(pde.fieldName);
        }
    }
    // model.result().export("data1").set("unit", units.toArray(new String[0]));
    buffer.append("model.result().export('data1').set('unit'," + tojs(units.toArray(new String[0])) + ");" + "\n");
    // model.result().export("data1").set("descr", descriptions.toArray(new String[0]));
    buffer.append("model.result().export('data1').set('descr'," + tojs(descriptions.toArray(new String[0])) + ");" + "\n");
    // model.result().export("data1").set("expr", expressions.toArray(new String[0]));
    buffer.append("model.result().export('data1').set('expr'," + tojs(expressions.toArray(new String[0])) + ");" + "\n");
    System.out.println("about to run() report");
    // model.result().export("data1").run();
    buffer.append("model.result().export('data1').run();" + "\n");
    // long t_endexport_ms = System.currentTimeMillis();
    buffer.append("var t_endexport_ms = Date.now();" + "\n");
    // System.out.println("export finshed in " + ((t_endexport_ms - t_endsolver_ms) / 1000.0) + " seconds");
    buffer.append("print('export finshed in ' + ((t_endexport_ms - t_endsolver_ms) / 1000.0) + ' seconds');" + "\n");
    // model.save(javaFile.getAbsolutePath(), "java");
    String javaFileName = StringEscapeUtils.escapeEcmaScript(javaFile.getAbsolutePath());
    buffer.append("model.save('" + javaFileName + "','java');" + "\n");
    // model.save(mphFile.getAbsolutePath());
    String mphFileName = StringEscapeUtils.escapeEcmaScript(mphFile.getAbsolutePath());
    buffer.append("model.save('" + mphFileName + "');" + "\n");
    // ModelUtil.disconnect();
    buffer.append("com.comsol.model.util.ModelUtil.disconnect();" + "\n");
    String script = buffer.toString();
    System.out.println(script);
    FileUtils.write(jsFile, script);
    System.out.println("script written to " + jsFile.getAbsolutePath() + " prior to execution");
    // RhinoScriptEngineFactory factory = new RhinoScriptEngineFactory();
    // factory.
    // ContextFactory factory = new ContextFactory();
    // Main dbg = new Main("script");
    // dbg.attachTo(factory);
    // Context context = factory.
    // .
    // ScriptContext scriptContext = new ScriptContext
    engine.eval(script);
}
Also used : VCCDifference(org.vcell.solver.comsol.model.VCCGeomFeature.VCCDifference) VCCGeomSequence(org.vcell.solver.comsol.model.VCCGeomSequence) VCCPhysicsFeature(org.vcell.solver.comsol.model.VCCPhysicsFeature) VCCMove(org.vcell.solver.comsol.model.VCCGeomFeature.VCCMove) VCCIntersectionSelection(org.vcell.solver.comsol.model.VCCGeomFeature.VCCIntersectionSelection) VCCBlock(org.vcell.solver.comsol.model.VCCGeomFeature.VCCBlock) VCCConvectionDiffusionEquation(org.vcell.solver.comsol.model.VCCConvectionDiffusionEquation) ArrayList(java.util.ArrayList) VCCSphere(org.vcell.solver.comsol.model.VCCGeomFeature.VCCSphere) VCCModelNode(org.vcell.solver.comsol.model.VCCModelNode) VCCSquare(org.vcell.solver.comsol.model.VCCGeomFeature.VCCSquare) VCCFluxBoundary(org.vcell.solver.comsol.model.VCCPhysicsFeature.VCCFluxBoundary) VCCTransientStudyFeature(org.vcell.solver.comsol.model.VCCTransientStudyFeature) VCCScale(org.vcell.solver.comsol.model.VCCGeomFeature.VCCScale) VCCMeshSequence(org.vcell.solver.comsol.model.VCCMeshSequence) VCCStudyFeature(org.vcell.solver.comsol.model.VCCStudyFeature) VCCCircle(org.vcell.solver.comsol.model.VCCGeomFeature.VCCCircle) VCCPhysics(org.vcell.solver.comsol.model.VCCPhysics) VCCGeomFeature(org.vcell.solver.comsol.model.VCCGeomFeature)

Example 3 with VCCDifference

use of org.vcell.solver.comsol.model.VCCGeomFeature.VCCDifference in project vcell by virtualcell.

the class ComsolModelBuilder method getVCCModel.

public static VCCModel getVCCModel(SimulationJob vcellSimJob) throws ExpressionException {
    MathDescription vcellMathDesc = vcellSimJob.getSimulation().getMathDescription();
    Geometry vcellGeometry = vcellMathDesc.getGeometry();
    GeometrySpec vcellGeometrySpec = vcellGeometry.getGeometrySpec();
    int vcellDim = vcellGeometrySpec.getDimension();
    VCCModel model = new VCCModel("Model", vcellDim);
    model.modelpath = "D:\\Developer\\eclipse\\workspace_refactor\\comsol_java\\src";
    model.comments = "Untitled\n\n";
    VCCModelNode comp1 = new VCCModelNode("comp1");
    model.modelnodes.add(comp1);
    // if (vcellDim != 2){
    // throw new RuntimeException("expecting 2D simulation");
    // }
    // 
    // assume initial geometry is circle centered at 0.5, 0.5 of radius 0.3
    // 
    // String comsolOutsideDomainName = "dif1";
    // String comsolInsideDomainName = "c1";
    VCCGeomSequence geom1 = new VCCGeomSequence("geom1", vcellDim);
    model.geometrysequences.add(geom1);
    VCCMeshSequence mesh1 = new VCCMeshSequence("mesh1", geom1);
    model.meshes.add(mesh1);
    VCCStudy std1 = new VCCStudy("std1");
    model.study = std1;
    TimeBounds timeBounds = vcellSimJob.getSimulation().getSolverTaskDescription().getTimeBounds();
    TimeStep timeStep = vcellSimJob.getSimulation().getSolverTaskDescription().getTimeStep();
    String beginTime = Double.toString(timeBounds.getStartingTime());
    String endTime = Double.toString(timeBounds.getEndingTime());
    String step = Double.toString(timeStep.getDefaultTimeStep());
    VCCStudyFeature time = new VCCTransientStudyFeature("time", beginTime, step, endTime);
    std1.features.add(time);
    if (vcellGeometrySpec.getImage() != null) {
        throw new RuntimeException("image-based geometries not yet supported by VCell's COMSOL model builder");
    }
    if (vcellGeometrySpec.getNumSubVolumes() == 0) {
        throw new RuntimeException("no subvolumes defined in geometry");
    }
    if (vcellGeometrySpec.getNumAnalyticOrCSGSubVolumes() != vcellGeometrySpec.getNumSubVolumes()) {
        throw new RuntimeException("only analytic and CSG subvolumes currently supported by VCell's COMSOL model builder");
    }
    // 
    // add geometry for all subvolumes
    // 
    HashMap<String, VCCGeomFeature> subvolumeNameFeatureMap = new HashMap<String, VCCGeomFeature>();
    SubVolume[] subVolumes = vcellGeometrySpec.getSubVolumes();
    for (int i = 0; i < subVolumes.length; i++) {
        SubVolume subvolume = subVolumes[i];
        if (subvolume instanceof CSGObject) {
            CSGObject vcellCSGObject = (CSGObject) subvolume;
            CSGNode vcellCSGNode = vcellCSGObject.getRoot();
            ArrayList<VCCGeomFeature> geomFeatureList = new ArrayList<VCCGeomFeature>();
            VCCGeomFeature feature = csgVisitor(vcellCSGNode, geomFeatureList, subvolume.getName());
            geom1.geomfeatures.addAll(geomFeatureList);
            if (i == 0) {
                // first subvolume (on top in ordinals) doesn't need any differencing
                subvolumeNameFeatureMap.put(subvolume.getName(), feature);
            } else {
                // have to subtract union of prior subvolumes
                ArrayList<VCCGeomFeature> priorFeatures = new ArrayList<VCCGeomFeature>();
                for (int j = 0; j < i; j++) {
                    CSGObject priorCSGObject = (CSGObject) subVolumes[j];
                    CSGNode priorCSGNode = priorCSGObject.getRoot();
                    geomFeatureList.clear();
                    VCCGeomFeature priorFeature = csgVisitor(priorCSGNode, geomFeatureList, subvolume.getName());
                    priorFeatures.add(priorFeature);
                    geom1.geomfeatures.addAll(geomFeatureList);
                }
                VCCDifference diff = new VCCDifference("diff" + subvolume.getName(), Keep.off);
                diff.input.add(feature);
                diff.input2.addAll(priorFeatures);
                geom1.geomfeatures.add(diff);
                subvolumeNameFeatureMap.put(subvolume.getName(), diff);
            }
        } else {
            throw new RuntimeException("only CSG subvolumes currently supported by VCell's COMSOL model builder");
        }
    }
    // 
    // add geometry for all surfaceClasses
    // 
    HashMap<String, VCCGeomFeature> surfaceclassNameFeatureMap = new HashMap<String, VCCGeomFeature>();
    SurfaceClass[] surfaceClasses = vcellGeometry.getGeometrySurfaceDescription().getSurfaceClasses();
    for (int i = 0; i < surfaceClasses.length; i++) {
        SurfaceClass surfaceClass = surfaceClasses[i];
        Set<SubVolume> adjacentSubvolumes = surfaceClass.getAdjacentSubvolumes();
        if (adjacentSubvolumes.size() != 2) {
            throw new RuntimeException("expecting two adjacent subvolumes for surface " + surfaceClass.getName() + " in COMSOL model builder");
        }
        // find adjacent Geometry Features (for subvolumes)
        Iterator<SubVolume> svIter = adjacentSubvolumes.iterator();
        SubVolume subvolume0 = svIter.next();
        SubVolume subvolume1 = svIter.next();
        ArrayList<VCCGeomFeature> adjacentFeatures = new ArrayList<VCCGeomFeature>();
        adjacentFeatures.add(subvolumeNameFeatureMap.get(subvolume0.getName()));
        adjacentFeatures.add(subvolumeNameFeatureMap.get(subvolume1.getName()));
        String name = "inter_" + subvolume0.getName() + "_" + subvolume1.getName();
        // surfaces are dimension N-1
        int entitydim = vcellDim - 1;
        VCCIntersectionSelection intersect_subvolumes = new VCCIntersectionSelection(name, entitydim);
        intersect_subvolumes.input.addAll(adjacentFeatures);
        geom1.geomfeatures.add(intersect_subvolumes);
        surfaceclassNameFeatureMap.put(surfaceClass.getName(), intersect_subvolumes);
    }
    SimulationSymbolTable symbolTable = new SimulationSymbolTable(vcellSimJob.getSimulation(), vcellSimJob.getJobIndex());
    // 
    for (SubDomain subDomain : Collections.list(vcellMathDesc.getSubDomains())) {
        for (Equation equ : subDomain.getEquationCollection()) {
            if (equ instanceof PdeEquation || equ instanceof OdeEquation) {
                VCCGeomFeature geomFeature = null;
                final int dim;
                if (subDomain instanceof CompartmentSubDomain) {
                    geomFeature = subvolumeNameFeatureMap.get(subDomain.getName());
                    dim = vcellDim;
                } else if (subDomain instanceof MembraneSubDomain) {
                    geomFeature = surfaceclassNameFeatureMap.get(subDomain.getName());
                    dim = vcellDim - 1;
                } else {
                    throw new RuntimeException("subdomains of type '" + subDomain.getClass().getSimpleName() + "' not yet supported in COMSOL model builder");
                }
                if (geomFeature == null) {
                    throw new RuntimeException("cannot find COMSOL geometry feature named " + subDomain.getName() + " in COMSOL model builder");
                }
                VCCConvectionDiffusionEquation cdeq = new VCCConvectionDiffusionEquation("cdeq_" + equ.getVariable().getName(), geom1, geomFeature, dim);
                cdeq.fieldName = equ.getVariable().getName();
                cdeq.initial = MathUtilities.substituteModelParameters(equ.getInitialExpression(), symbolTable).flatten().infix();
                cdeq.sourceTerm_f = MathUtilities.substituteModelParameters(equ.getRateExpression(), symbolTable).flatten().infix();
                if (equ instanceof PdeEquation) {
                    PdeEquation pde = (PdeEquation) equ;
                    cdeq.diffTerm_c = MathUtilities.substituteModelParameters(pde.getDiffusionExpression(), symbolTable).flatten().infix();
                    if (subDomain instanceof CompartmentSubDomain) {
                        CompartmentSubDomain compartmentSubdomain = (CompartmentSubDomain) subDomain;
                        ArrayList<String> be = new ArrayList<String>();
                        if (pde.getVelocityX() != null) {
                            be.add(MathUtilities.substituteModelParameters(pde.getVelocityX(), symbolTable).flatten().infix());
                        } else {
                            be.add("0");
                        }
                        if (vcellDim >= 2) {
                            if (pde.getVelocityY() != null) {
                                be.add(MathUtilities.substituteModelParameters(pde.getVelocityY(), symbolTable).flatten().infix());
                            } else {
                                be.add("0");
                            }
                        }
                        if (vcellDim == 3) {
                            if (pde.getVelocityY() != null) {
                                be.add(MathUtilities.substituteModelParameters(pde.getVelocityZ(), symbolTable).flatten().infix());
                            } else {
                                be.add("0");
                            }
                        }
                        cdeq.advection_be = be.toArray(new String[vcellDim]);
                        // 
                        // look for membrane boundary conditions for this variable
                        // 
                        MembraneSubDomain[] membraneSubdomains = vcellMathDesc.getMembraneSubDomains(compartmentSubdomain);
                        for (MembraneSubDomain membraneSubdomain : membraneSubdomains) {
                            JumpCondition jumpCondition = membraneSubdomain.getJumpCondition((VolVariable) pde.getVariable());
                            if (jumpCondition != null) {
                                Expression fluxExpr = null;
                                if (membraneSubdomain.getInsideCompartment() == compartmentSubdomain) {
                                    fluxExpr = jumpCondition.getInFluxExpression();
                                } else if (membraneSubdomain.getOutsideCompartment() == compartmentSubdomain) {
                                    fluxExpr = jumpCondition.getOutFluxExpression();
                                }
                                String name = equ.getVariable().getName() + "_flux_" + membraneSubdomain.getName();
                                VCCGeomFeature selection = surfaceclassNameFeatureMap.get(membraneSubdomain.getName());
                                VCCFluxBoundary fluxBoundary = new VCCFluxBoundary(name, selection, vcellDim - 1);
                                fluxBoundary.flux_g = MathUtilities.substituteModelParameters(fluxExpr, symbolTable).flatten().infix();
                                cdeq.features.add(fluxBoundary);
                            }
                        }
                    }
                }
                model.physics.add(cdeq);
            }
        }
    }
    // 
    return model;
}
Also used : JumpCondition(cbit.vcell.math.JumpCondition) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) MathDescription(cbit.vcell.math.MathDescription) HashMap(java.util.HashMap) SurfaceClass(cbit.vcell.geometry.SurfaceClass) VCCConvectionDiffusionEquation(org.vcell.solver.comsol.model.VCCConvectionDiffusionEquation) CSGNode(cbit.vcell.geometry.CSGNode) ArrayList(java.util.ArrayList) VCCFluxBoundary(org.vcell.solver.comsol.model.VCCPhysicsFeature.VCCFluxBoundary) GeometrySpec(cbit.vcell.geometry.GeometrySpec) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) PdeEquation(cbit.vcell.math.PdeEquation) TimeBounds(cbit.vcell.solver.TimeBounds) TimeStep(cbit.vcell.solver.TimeStep) SubVolume(cbit.vcell.geometry.SubVolume) VCCStudyFeature(org.vcell.solver.comsol.model.VCCStudyFeature) CSGObject(cbit.vcell.geometry.CSGObject) VCCGeomFeature(org.vcell.solver.comsol.model.VCCGeomFeature) VCCDifference(org.vcell.solver.comsol.model.VCCGeomFeature.VCCDifference) VCCGeomSequence(org.vcell.solver.comsol.model.VCCGeomSequence) VCCIntersectionSelection(org.vcell.solver.comsol.model.VCCGeomFeature.VCCIntersectionSelection) SimulationSymbolTable(cbit.vcell.solver.SimulationSymbolTable) VCCConvectionDiffusionEquation(org.vcell.solver.comsol.model.VCCConvectionDiffusionEquation) OdeEquation(cbit.vcell.math.OdeEquation) PdeEquation(cbit.vcell.math.PdeEquation) Equation(cbit.vcell.math.Equation) VCCModelNode(org.vcell.solver.comsol.model.VCCModelNode) VCCTransientStudyFeature(org.vcell.solver.comsol.model.VCCTransientStudyFeature) VCCMeshSequence(org.vcell.solver.comsol.model.VCCMeshSequence) Geometry(cbit.vcell.geometry.Geometry) VCCStudy(org.vcell.solver.comsol.model.VCCStudy) OdeEquation(cbit.vcell.math.OdeEquation) Expression(cbit.vcell.parser.Expression) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) VCCModel(org.vcell.solver.comsol.model.VCCModel)

Aggregations

ArrayList (java.util.ArrayList)3 VCCGeomFeature (org.vcell.solver.comsol.model.VCCGeomFeature)3 VCCDifference (org.vcell.solver.comsol.model.VCCGeomFeature.VCCDifference)3 CSGNode (cbit.vcell.geometry.CSGNode)2 VCCConvectionDiffusionEquation (org.vcell.solver.comsol.model.VCCConvectionDiffusionEquation)2 VCCBlock (org.vcell.solver.comsol.model.VCCGeomFeature.VCCBlock)2 VCCIntersectionSelection (org.vcell.solver.comsol.model.VCCGeomFeature.VCCIntersectionSelection)2 VCCMove (org.vcell.solver.comsol.model.VCCGeomFeature.VCCMove)2 VCCScale (org.vcell.solver.comsol.model.VCCGeomFeature.VCCScale)2 VCCSphere (org.vcell.solver.comsol.model.VCCGeomFeature.VCCSphere)2 VCCGeomSequence (org.vcell.solver.comsol.model.VCCGeomSequence)2 VCCMeshSequence (org.vcell.solver.comsol.model.VCCMeshSequence)2 VCCModelNode (org.vcell.solver.comsol.model.VCCModelNode)2 VCCFluxBoundary (org.vcell.solver.comsol.model.VCCPhysicsFeature.VCCFluxBoundary)2 VCCStudyFeature (org.vcell.solver.comsol.model.VCCStudyFeature)2 VCCTransientStudyFeature (org.vcell.solver.comsol.model.VCCTransientStudyFeature)2 CSGHomogeneousTransformation (cbit.vcell.geometry.CSGHomogeneousTransformation)1 CSGObject (cbit.vcell.geometry.CSGObject)1 CSGPrimitive (cbit.vcell.geometry.CSGPrimitive)1 CSGRotation (cbit.vcell.geometry.CSGRotation)1