/*
 * cohesive_zone_debonding.java
 */

import com.comsol.model.*;
import com.comsol.model.util.*;

/** Model exported on May 13 2026, 11:26 by COMSOL 6.4.0.419. */
public class cohesive_zone_debonding {

  public static Model run() {
    Model model = ModelUtil.create("Model");

//    From the File menu, choose New.
//    In the New window, click Model Wizard.
//    In the Model Wizard window, click 3D.
//    In the Select Physics tree, select Structural Mechanics > Solid Mechanics (solid).
//    Click Add.
//    Click Study.
//    In the Select Study tree, select General Studies > Stationary.
//    Click Done.

    model.component().create("comp1", true);

    model.component("comp1").geom().create("geom1", 3);
    model.component("comp1").geom("geom1").geomRep("comsol");

    model.component("comp1").mesh().create("mesh1");
    model.component("comp1").mesh("mesh1").contribute("geom/detail", true);

    model.component("comp1").physics().create("solid", "SolidMechanics", "geom1");

    model.study().create("std1");
    model.study("std1").create("stat", "Stationary");

//    Load all model parameters from a file containing parameters for the geometry, material properties and boundary conditions.
//    In the Model Builder window, under Global Definitions, click Parameters 1.
//    In the Settings window for Parameters, locate the Parameters section.
//    Click Load from File.
//    Browse to the model's Application Library folder and double-click the file cohesive_zone_debonding_parameters.txt.
//    To import content from file, use:
//    model.param().loadFile("FILENAME");
    model.param().set("lb", "102[mm]", "Length");
    model.param().set("wb", "25.4[mm]", "Width");
    model.param().set("hb", "2*1.56[mm]", "Thickness");
    model.param().set("cl", "34.1[mm]", "Initial crack length");
    model.param().set("pn", "1e6[N/mm^3]", "Penalty stiffness");
    model.param().set("sigmat", "80[MPa]", "Normal tensile strength");
    model.param().set("sigmas", "100[MPa]", "Shear strength");
    model.param().set("Gct", "0.969[kJ/m^2]", "Critical energy release rate, tension");
    model.param().set("Gcs", "1.719[kJ/m^2]", "Critical energy release rate, shear");
    model.param().set("alpha", "2.284", "Exponent of the Benzeggagh and Kenane (B-K) criterion");
    model.param().set("disp", "0", "Displacement parameter");
    model.param().set("mm", "0.5", "Mode mixity ratio");
    model.param().set("ll", "lb/2*(0.5*sqrt(3*(1-mm)/mm)+1)/(3-0.5*sqrt(3*(1-mm)/mm))", "Lever length");
    model.param()
         .set("lr", "8*((6*mm+sqrt(3*mm*(1-mm)))/(3+9*mm+8*sqrt(3*mm*(1-mm))))", "Load ratio middle/cracked edge");

//    In the Model Builder window, under Component 1 (comp1), click Geometry 1.
//    In the Settings window for Geometry, locate the Units section.
//    From the Length unit list, select mm.

    model.component("comp1").geom("geom1").lengthUnit("mm");

//    Now draw the model geometry with two layered blocks that are overlapping.
//    In the Geometry toolbar, click Block.

    model.component("comp1").geom("geom1").create("blk1", "Block");

//    In the Settings window for Block, locate the Size and Shape section.
//    In the Width text field, type lb.

    model.component("comp1").geom("geom1").feature("blk1").set("size", new String[]{"lb", "1", "1"});

//    In the Depth text field, type wb/2.

    model.component("comp1").geom("geom1").feature("blk1").set("size", new String[]{"lb", "wb/2", "1"});

//    In the Height text field, type hb.

    model.component("comp1").geom("geom1").feature("blk1").set("size", new String[]{"lb", "wb/2", "hb"});

//    Click to expand the Layers section.
//    In the table, enter the following settings:

    model.component("comp1").geom("geom1").feature("blk1").setIndex("layer", "hb/2", 0);

//    In the Geometry toolbar, click Block.

    model.component("comp1").geom("geom1").run("blk1");
    model.component("comp1").geom("geom1").create("blk2", "Block");

//    In the Settings window for Block, locate the Size and Shape section.
//    In the Width text field, type lb/2-cl.

    model.component("comp1").geom("geom1").feature("blk2").set("size", new String[]{"lb/2-cl", "1", "1"});

//    In the Depth text field, type wb/2.

    model.component("comp1").geom("geom1").feature("blk2").set("size", new String[]{"lb/2-cl", "wb/2", "1"});

//    In the Height text field, type hb.

    model.component("comp1").geom("geom1").feature("blk2").setIndex("size", "hb", 2);

//    Locate the Position section.
//    In the x text field, type cl.

    model.component("comp1").geom("geom1").feature("blk2").set("pos", new String[]{"cl", "0", "0"});

//    Locate the Layers section.
//    In the table, enter the following settings:

    model.component("comp1").geom("geom1").feature("blk2").setIndex("layer", "hb/2", 0);

//    In the Geometry toolbar, click Build All.

    model.component("comp1").geom("geom1").run("fin");

//    In the Model Builder window, expand the Component 1 (comp1) > Definitions node, then click Boundary System 1 (sys1).
//    In the Settings window for Boundary System, locate the Settings section.
//    From the Frame list, select Reference configuration.

    model.component("comp1").coordSystem("sys1").set("frametype", "material");

//    In the Definitions toolbar, click Local Variables.

    model.component("comp1").variable().create("var1");

//    Load variables for the load point from files.
//    In the Settings window for Variables, type Load Point Variables in the Label text field.

    model.component("comp1").variable("var1").label("Load Point Variables");

//    Locate the Variables section.
//    Click Load from File.
//    Browse to the model's Application Library folder and double-click the file cohesive_zone_debonding_load_point_variables.txt.
//    To import content from file, use:
//    model.component("comp1").variable("var1").loadFile("FILENAME");
    model.component("comp1").variable("var1").set("u_Ie", "intop1(w)", "Displacement at edge");
    model.component("comp1").variable("var1").set("w_c", "intop2(w)", "Displacement at center");
    model.component("comp1").variable("var1")
         .set("u_lp", "(3*ll-lb/2)/4/(lb/2)*u_Ie+((ll+lb/2)/(lb/2))*(-w_c+u_Ie/4)", "Load point displacement");
    model.component("comp1").variable("var1").set("F_lp", "force*lb/2/ll", "Load point force");

//    In the Definitions toolbar, click Nonlocal Couplings and choose Integration.

    model.component("comp1").cpl().create("intop1", "Integration");
    model.component("comp1").cpl("intop1").set("axisym", true);

//    The following nonlocal integration couplings make values of the selected points globally available.
//    In the Settings window for Integration, type Integration Edge in the Label text field.

    model.component("comp1").cpl("intop1").label("Integration Edge");

//    Locate the Source Selection section.
//    From the Geometric entity level list, select Point.

    model.component("comp1").cpl("intop1").selection().geom("geom1", 0);

//    Select Point 3.

    model.component("comp1").cpl("intop1").selection().set(3);

//    Locate the Advanced section.
//    From the Frame list, select Material (X, Y, Z).

    model.component("comp1").cpl("intop1").set("frame", "material");

//    In the Definitions toolbar, click Nonlocal Couplings and choose Integration.

    model.component("comp1").cpl().create("intop2", "Integration");
    model.component("comp1").cpl("intop2").set("axisym", true);

//    In the Settings window for Integration, type Integration Center in the Label text field.

    model.component("comp1").cpl("intop2").label("Integration Center");

//    Locate the Source Selection section.
//    From the Geometric entity level list, select Point.

    model.component("comp1").cpl("intop2").selection().geom("geom1", 0);

//    Select Point 15.

    model.component("comp1").cpl("intop2").selection().set(15);

//    Locate the Advanced section.
//    From the Frame list, select Material (X, Y, Z).

    model.component("comp1").cpl("intop2").set("frame", "material");

//    COMSOL Multiphysics is equipped with built-in material properties for a number of lamina materials. Select the needed materials from the <l>Composites</l> material folder in the built-in material library.
//    In the Materials toolbar, click Add Material to open the Add Material window.
//    In the tree, select Composites > Laminae > Unidirectional fiber lamina: AS4/APC2 carbon/PEEK thermoplastic [fiber volume fraction 50%].
//    Right-click and choose Add to Global Materials.

    model.material().create("mat1", "Common", "");
    model.material("mat1").propertyGroup()
         .create("TransverseIsotropic", "TransverseIsotropic", "Transversely isotropic");
    model.material("mat1")
         .label("Unidirectional fiber lamina: AS4/APC2 carbon/PEEK thermoplastic [fiber volume fraction 50%]");
    model.material("mat1").propertyGroup("def").set("density", "1570[kg/m^3]");
    model.material("mat1").propertyGroup("def")
         .setPropertyInfo("density", "Reference: P.P. Camanho, C.G. Davila, and M.F. De Moura, Numerical Simulation of Mixed Mode Progressive Delamination in Composite Materials, Journal of composite materials, vol. 37, no. 16, 2003.");
    model.material("mat1").propertyGroup("def")
         .set("thermalexpansioncoefficient", new String[]{"-0.2E-6[1/K]", "0", "0", "0", "24E-6[1/K]", "0", "0", "0", "24E-6[1/K]"});
    model.material("mat1").propertyGroup("TransverseIsotropic").set("Evect", new String[]{"122.7[GPa]", "10.1[GPa]"});
    model.material("mat1").propertyGroup("TransverseIsotropic")
         .setPropertyInfo("Evect", "Reference: P.P. Camanho, C.G. Davila, and M.F. De Moura, Numerical Simulation of Mixed Mode Progressive Delamination in Composite Materials, Journal of composite materials, vol. 37, no. 16, 2003.");
    model.material("mat1").propertyGroup("TransverseIsotropic").set("nuvect", new String[]{"0.25", "0.45"});
    model.material("mat1").propertyGroup("TransverseIsotropic")
         .setPropertyInfo("nuvect", "Reference: P.P. Camanho, C.G. Davila, and M.F. De Moura, Numerical Simulation of Mixed Mode Progressive Delamination in Composite Materials, Journal of composite materials, vol. 37, no. 16, 2003.");
    model.material("mat1").propertyGroup("TransverseIsotropic").set("Gvect1", "5.5[GPa]");
    model.material("mat1").propertyGroup("TransverseIsotropic")
         .setPropertyInfo("Gvect1", "Reference: P.P. Camanho, C.G. Davila, and M.F. De Moura, Numerical Simulation of Mixed Mode Progressive Delamination in Composite Materials, Journal of composite materials, vol. 37, no. 16, 2003.");

//    In the Materials toolbar, click Add Material to close the Add Material window.
//    In the Model Builder window, right-click Component 1 (comp1) > Materials and choose More Materials > Material Link.

    model.component("comp1").material().create("matlnk1", "Link");

//    In the Settings window for Linear Elastic Material, locate the Linear Elastic Material section.
//    From the Material symmetry list, select Orthotropic.

    model.component("comp1").physics("solid").feature("lemm1").set("SolidModel", "Orthotropic");

//    Add the CZM using a <l>Interior Contact</l> node.
//    In the Physics toolbar, click Boundaries and choose Interior Contact.

    model.component("comp1").physics("solid").create("icnt1", "InteriorContact", 2);

//    Select Boundaries 15, 24.

    model.component("comp1").physics("solid").feature("icnt1").selection().set(15, 24);

//    In the Settings window for Interior Contact, locate the Contact Pressure Penalty Factor section.
//    From the Penalty factor control list, select User defined.

    model.component("comp1").physics("solid").feature("icnt1").set("penaltyCtrlPenalty", "userDefined");

//    In the \[p_{\textrm{n}}\] text field, type pn.

    model.component("comp1").physics("solid").feature("icnt1").set("pn_penalty", "pn");

//    In the Physics toolbar, click Attributes and choose Adhesion.

    model.component("comp1").physics("solid").feature("icnt1").create("adh1", "Adhesion", 2);

//    In the Settings window for Adhesion, locate the Adhesive Activation section.
//    From the Activation criterion list, select Always active.

    model.component("comp1").physics("solid").feature("icnt1").feature("adh1").set("ActivationCriterion", "Always");

//    Locate the Adhesive Stiffness section.
//    In the \[n_{\tau}\] text field, type 1.

    model.component("comp1").physics("solid").feature("icnt1").feature("adh1").set("ntau", 1);

//    In the Physics toolbar, click Attributes and choose Decohesion.

    model.component("comp1").physics("solid").feature("icnt1").feature("adh1").create("dch1", "Decohesion", 2);

//    In the Settings window for Decohesion, locate the Decohesion section.
//    In the \[\sigma_{\textrm{ts}}\] text field, type sigmat.

    model.component("comp1").physics("solid").feature("icnt1").feature("adh1").feature("dch1")
         .set("sigmat", "sigmat");

//    In the \[\sigma_{\textrm{ss}}\] text field, type sigmas.

    model.component("comp1").physics("solid").feature("icnt1").feature("adh1").feature("dch1")
         .set("sigmas", "sigmas");

//    In the \[G_{\textrm{ct}}\] text field, type Gct.

    model.component("comp1").physics("solid").feature("icnt1").feature("adh1").feature("dch1").set("Gct", "Gct");

//    In the \[G_{\textrm{cs}}\] text field, type Gcs.

    model.component("comp1").physics("solid").feature("icnt1").feature("adh1").feature("dch1").set("Gcs", "Gcs");

//    From the Mixed mode criterion list, select Benzeggagh–Kenane.

    model.component("comp1").physics("solid").feature("icnt1").feature("adh1").feature("dch1")
         .set("FailureCriterion", "BK");

//    In the \[\alpha\] text field, type alpha.

    model.component("comp1").physics("solid").feature("icnt1").feature("adh1").feature("dch1").set("alpha", "alpha");

//    Click the Show More Options button in the Model Builder toolbar.
//    In the Show More Options dialog, in the tree, select the checkbox for the Physics > Equation Contributions node.
//    Click OK.
//    This is to make <l>Auxiliary Slit</l> and <l>Global Equations</l> accessible.
//    In the Physics toolbar, click Boundaries and choose Auxiliary Slit.

    model.component("comp1").physics("solid").create("asl1", "AuxiliarySlit", 2);

//    Select Boundary 6.

    model.component("comp1").physics("solid").feature("asl1").selection().set(6);

//    In the Physics toolbar, click Boundaries and choose Symmetry.

    model.component("comp1").physics("solid").create("sym1", "SymmetrySolid", 2);

//    Select Boundaries 2, 5, 11, 14, 20, 23.

    model.component("comp1").physics("solid").feature("sym1").selection().set(2, 5, 11, 14, 20, 23);

//    In the Physics toolbar, click Edges and choose Edge Load.

    model.component("comp1").physics("solid").create("el1", "EdgeLoad", 1);

//    In the Settings window for Edge Load, type Load on Cracked Edge (Fe) in the Label text field.

    model.component("comp1").physics("solid").feature("el1").label("Load on Cracked Edge (Fe)");

//    Select Edge 7.

    model.component("comp1").physics("solid").feature("el1").selection().set(7);

//    Locate the Force section.
//    From the Load type list, select Total force.

    model.component("comp1").physics("solid").feature("el1").set("forceType", "TotalForce");

//    Specify the \[\mathbf{F}_{\mathrm{tot}}\] vector as

    model.component("comp1").physics("solid").feature("el1").set("force", new String[]{"0", "0", "force"});

//    In the Physics toolbar, click Edges and choose Edge Load.

    model.component("comp1").physics("solid").create("el2", "EdgeLoad", 1);

//    In the Settings window for Edge Load, type Load on Middle Edge (Fm) in the Label text field.

    model.component("comp1").physics("solid").feature("el2").label("Load on Middle Edge (Fm)");

//    Select Edge 33.

    model.component("comp1").physics("solid").feature("el2").selection().set(33);

//    Locate the Force section.
//    From the Load type list, select Total force.

    model.component("comp1").physics("solid").feature("el2").set("forceType", "TotalForce");

//    Specify the \[\mathbf{F}_{\mathrm{tot}}\] vector as

    model.component("comp1").physics("solid").feature("el2").set("force", new String[]{"0", "0", "-lr*force"});

//    In the Physics toolbar, click Edges and choose Prescribed Displacement.

    model.component("comp1").physics("solid").create("disp1", "Displacement1", 1);

//    Select Edges 2, 41.

    model.component("comp1").physics("solid").feature("disp1").selection().set(2, 41);

//    In the Settings window for Prescribed Displacement, locate the Prescribed Displacement section.
//    From the Displacement in z direction list, select Prescribed.

    model.component("comp1").physics("solid").feature("disp1").setIndex("Direction", "prescribed", 2);

//    In the Physics toolbar, click Points and choose Prescribed Displacement.

    model.component("comp1").physics("solid").create("disp2", "Displacement0", 0);

//    Select Point 1.

    model.component("comp1").physics("solid").feature("disp2").selection().set(1);

//    In the Settings window for Prescribed Displacement, locate the Prescribed Displacement section.
//    From the Displacement in x direction list, select Prescribed.

    model.component("comp1").physics("solid").feature("disp2").setIndex("Direction", "prescribed", 0);

//    Add a global equation to control the applied load with a monotonically increasing parameter.
//    In the Physics toolbar, click Global and choose Global Equations.

    model.component("comp1").physics("solid").create("ge1", "GlobalEquations", -1);

//    In the Settings window for Global Equations, locate the Global Equations section.
//    In the table, enter the following settings:

    model.component("comp1").physics("solid").feature("ge1").setIndex("name", "force", 0, 0);
    model.component("comp1").physics("solid").feature("ge1").setIndex("equation", "disp-u_Ie", 0, 0);

//    Locate the Units section.
//    Click Select Dependent Variable Quantity.
//    In the Physical Quantity dialog, type force in the text field.
//    In the tree, select General > Force (N).
//    Click OK.

    model.component("comp1").physics("solid").feature("ge1").set("DependentVariableQuantity", "force");

//    In the Settings window for Global Equations, locate the Units section.
//    Click Select Source Term Quantity.
//    In the Physical Quantity dialog, type Displacement in the text field.
//    In the tree, select General > Displacement (m).
//    Click OK.

    model.component("comp1").physics("solid").feature("ge1").set("SourceTermQuantity", "displacement");

//    In the Mesh toolbar, click More Generators and choose Mapped.

    model.component("comp1").mesh("mesh1").create("map1", "Map");

//    Select Boundaries 2, 5, 11, 14, 20, 23.

    model.component("comp1").mesh("mesh1").feature("map1").selection().set(2, 5, 11, 14, 20, 23);

//    In the Mesh toolbar, click Distribution.

    model.component("comp1").mesh("mesh1").feature("map1").create("dis1", "Distribution");

//    Select Edges 3, 6, 8.

    model.component("comp1").mesh("mesh1").feature("map1").feature("dis1").selection().set(3, 6, 8);

//    In the Settings window for Distribution, locate the Distribution section.
//    In the Number of elements text field, type 10.

    model.component("comp1").mesh("mesh1").feature("map1").feature("dis1").set("numelem", 10);

//    In the Mesh toolbar, click Distribution.

    model.component("comp1").mesh("mesh1").feature("map1").create("dis2", "Distribution");

//    Select Edges 16, 19, 21.

    model.component("comp1").mesh("mesh1").feature("map1").feature("dis2").selection().set(16, 19, 21);

//    In the Settings window for Distribution, locate the Distribution section.
//    In the Number of elements text field, type 50.

    model.component("comp1").mesh("mesh1").feature("map1").feature("dis2").set("numelem", 50);

//    In the Mesh toolbar, click Distribution.

    model.component("comp1").mesh("mesh1").feature("map1").create("dis3", "Distribution");

//    Select Edges 29, 32, 34.

    model.component("comp1").mesh("mesh1").feature("map1").feature("dis3").selection().set(29, 32, 34);

//    In the Settings window for Distribution, locate the Distribution section.
//    From the Distribution type list, select Predefined.

    model.component("comp1").mesh("mesh1").feature("map1").feature("dis3").set("type", "predefined");

//    In the Number of elements text field, type 20.

    model.component("comp1").mesh("mesh1").feature("map1").feature("dis3").set("elemcount", 20);

//    In the Element ratio text field, type 10.

    model.component("comp1").mesh("mesh1").feature("map1").feature("dis3").set("elemratio", 10);

//    In the Mesh toolbar, click Swept.

    model.component("comp1").mesh("mesh1").create("swe1", "Sweep");

//    In the Mesh toolbar, click Distribution.

    model.component("comp1").mesh("mesh1").feature("swe1").create("dis1", "Distribution");

//    In the Settings window for Distribution, locate the Distribution section.
//    In the Number of elements text field, type 3.

    model.component("comp1").mesh("mesh1").feature("swe1").feature("dis1").set("numelem", 3);

//    Click Build All.

    model.component("comp1").mesh("mesh1").run();

//    Configure the solver to enable tracking of the post peak behavior of the beam.
//    In the Model Builder window, under Study 1, click Step 1: Stationary.
//    In the Settings window for Stationary, locate the Study Settings section.
//    Select the Include geometric nonlinearity checkbox.

    model.study("std1").feature("stat").set("geometricNonlinearity", true);

//    From the Tolerance list, select User controlled.

    model.study("std1").feature("stat").set("usestol", true);

//    In the Relative tolerance text field, type 1e-4.

    model.study("std1").feature("stat").set("stol", "1e-4");

//    Click to expand the Study Extensions section.
//    Select the Auxiliary sweep checkbox.

    model.study("std1").feature("stat").set("useparam", true);

//    Click Add.

    model.study("std1").feature("stat").setIndex("pname", "lb", 0);
    model.study("std1").feature("stat").setIndex("plistarr", "", 0);
    model.study("std1").feature("stat").setIndex("punit", "m", 0);
    model.study("std1").feature("stat").setIndex("pname", "lb", 0);
    model.study("std1").feature("stat").setIndex("plistarr", "", 0);
    model.study("std1").feature("stat").setIndex("punit", "m", 0);

//    In the table, enter the following settings:

    model.study("std1").feature("stat").setIndex("pname", "disp", 0);
    model.study("std1").feature("stat").setIndex("plistarr", "range(0,2e-4,8e-3)", 0);

//    Click to expand the Results While Solving section.
//    Select the Plot checkbox.

    model.study("std1").feature("stat").set("plot", true);

//    In the Study toolbar, click Show Default Solver.

    model.study("std1").showAutoSequences("all");

//    In the Model Builder window, expand the Solution 1 (sol1) node.
//    In the Model Builder window, expand the Study 1 > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 node, then click Global Equations 1 (comp1.ODE1).
//    In the Settings window for State, locate the Scaling section.
//    From the Method list, select Manual.

    model.sol("sol1").feature("v1").feature("comp1_ODE1").set("scalemethod", "manual");

//    In the Scale text field, type 200.

    model.sol("sol1").feature("v1").feature("comp1_ODE1").set("scaleval", 200);

//    Use a linear predictor.
//    In the Model Builder window, expand the Study 1 > Solver Configurations > Solution 1 (sol1) > Stationary Solver 1 node, then click Parametric 1.
//    In the Settings window for Parametric, click to expand the Continuation section.
//    Select the Tuning of step size checkbox.

    model.sol("sol1").feature("s1").feature("p1").set("paramtuning", true);

//    In the Minimum step size text field, type 1e-6.

    model.sol("sol1").feature("s1").feature("p1").set("pminstep", "1e-6");

//    From the Predictor list, select Linear.

    model.sol("sol1").feature("s1").feature("p1").set("porder", "linear");

//    Switch to an undamped Newton method.
//    In the Model Builder window, under Study 1 > Solver Configurations > Solution 1 (sol1) > Stationary Solver 1, click Fully Coupled 1.
//    In the Settings window for Fully Coupled, click to expand the Method and Termination section.
//    From the Nonlinear method list, select Constant (Newton).

    model.sol("sol1").feature("s1").feature("fc1").set("dtech", "const");

//    In the Study toolbar, click Get Initial Value.

    model.study("std1").createAutoSequences("sol");
    model.study("std1").createAutoSequences("jobs");

    model.sol("sol1").runFromTo("st1", "v1");

    model.result().remove("pg1");

    model.study("std1").feature("stat").set("plotgroup", "Default");

    model.result().create("pg1", "PlotGroup3D");
    model.result("pg1").set("data", "dset1");
    model.result("pg1").setIndex("looplevel", 1, 0);
    model.result("pg1").label("Stress (solid)");
    model.result("pg1").set("frametype", "spatial");
    model.result("pg1").create("vol1", "Volume");
    model.result("pg1").feature("vol1").set("expr", new String[]{"solid.misesGp"});
    model.result("pg1").feature("vol1").set("threshold", "manual");
    model.result("pg1").feature("vol1").set("thresholdvalue", 0.2);
    model.result("pg1").feature("vol1").set("colortable", "Rainbow");
    model.result("pg1").feature("vol1").set("colortabletrans", "none");
    model.result("pg1").feature("vol1").set("colorscalemode", "linear");
    model.result("pg1").feature("vol1").set("resolution", "custom");
    model.result("pg1").feature("vol1").set("refine", 2);
    model.result("pg1").feature("vol1").set("colortable", "Prism");
    model.result("pg1").feature("vol1").create("def", "Deform");
    model.result("pg1").feature("vol1").feature("def").set("scaleactive", true);
    model.result("pg1").feature("vol1").feature("def").set("scale", "1");
    model.result("pg1").feature("vol1").feature("def").set("expr", new String[]{"u", "v", "w"});
    model.result("pg1").feature("vol1").feature("def").set("descr", "Displacement field");
    model.result("pg1").run();
    model.result("pg1").run();

//    In the Model Builder window, expand the Stress (solid) node, then click Volume 1.
//    In the Settings window for Volume, locate the Expression section.
//    From the Unit list, select MPa.

    model.result("pg1").feature("vol1").set("unit", "MPa");

//    The Internal Contact and Auxiliary Slit conditions will enforce different displacements on the top and the bottom sides of the same boundary. This will cause some plots with deformation to look strange. To remedy this, do the following:
//    Click to expand the Shrink Elements section.
//    In the Element scale factor text field, type 0.999.

    model.result("pg1").feature("vol1").set("elemscale", 0.999);
    model.result("pg1").run();

//    In the Model Builder window, click Stress (solid).
//    In the Settings window for 3D Plot Group, locate the Plot Settings section.
//    Clear the Plot dataset edges checkbox.

    model.result("pg1").set("edges", false);

//    In the Stress (solid) toolbar, click Plot.

    model.result("pg1").run();

//    In the Home toolbar, click Compute.

    model.study("std1").createAutoSequences("all");

    model.sol("sol1").runAll();

    model.result("pg1").run();

//    Now, add a plot showing the state of debonding at the laminate interface on the undeformed geometry.
//    In the Home toolbar, click Add Plot Group and choose 3D Plot Group.

    model.result().create("pg2", "PlotGroup3D");
    model.result("pg2").run();

//    In the Settings window for 3D Plot Group, type Interface Health in the Label text field.

    model.result("pg2").label("Interface Health");

//    Right-click Interface Health and choose Surface.

    model.result("pg2").create("surf1", "Surface");
    model.result("pg2").feature("surf1").set("evaluationsettings", "parent");

//    In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section.
//    From the menu, choose Component 1 (comp1) > Solid Mechanics > Contact > Adhesion > solid.bdmg - Damage - 1.

    model.result("pg2").feature("surf1").set("expr", "solid.bdmg");
    model.result("pg2").feature("surf1").set("descr", "Damage");

//    Locate the Coloring and Style section.
//    From the Color table list, select Traffic.

    model.result("pg2").feature("surf1").set("colortable", "Traffic");

//    This choice plots the debonded part in red while the healthy part remains green.
//    In the Interface Health toolbar, click Plot.

    model.result("pg2").run();

//    Since the changes in damage level are rather abrupt, a surface plot can be difficult to interpret. The actual evaluation of damage states during solution is done at the integration point level. Add a plot which displays the values in each Gauss point.
//    In the Results toolbar, click 3D Plot Group.

    model.result().create("pg3", "PlotGroup3D");
    model.result("pg3").run();

//    In the Settings window for 3D Plot Group, type Interface Health: Gauss Points in the Label text field.

    model.result("pg3").label("Interface Health: Gauss Points");

//    In the Interface Health: Gauss Points toolbar, click More Plots and choose Point.

    model.result("pg3").create("pt1", "Point");

//    In the Settings window for Point, locate the Expression section.
//    In the Expression text field, type solid.bdmg.

    model.result("pg3").feature("pt1").set("expr", "solid.bdmg");

//    Locate the Evaluation section.
//    From the Geometry level list, select Surface.

    model.result("pg3").feature("pt1").set("level", "surface");

//    From the Placement list, select Gauss points.

    model.result("pg3").feature("pt1").set("placement", "gausspoints");

//    In the Gauss-point order text field, type 4.

    model.result("pg3").feature("pt1").set("gporder", "4");

//    Locate the Coloring and Style section.
//    From the Color table list, select Traffic.

    model.result("pg3").feature("pt1").set("colortable", "Traffic");

//    In the Interface Health: Gauss Points toolbar, click Plot.

    model.result("pg3").run();

//    Select the Radius scale factor checkbox.

    model.result("pg3").feature("pt1").set("sphereradiusscaleactive", true);

//    In the associated text field, type 0.03.

    model.result("pg3").feature("pt1").set("sphereradiusscale", 0.03);

//    In the Interface Health: Gauss Points toolbar, click Plot.

    model.result("pg3").run();
    model.result("pg3").run();

//    In the Model Builder window, click Interface Health: Gauss Points.
//    In the Settings window for 3D Plot Group, locate the Plot Settings section.
//    Clear the Plot dataset edges checkbox.

    model.result("pg3").set("edges", false);

//    Right-click Interface Health: Gauss Points and choose Mesh.

    model.result("pg3").create("mesh1", "Mesh");
    model.result("pg3").feature("mesh1").set("colortable", "TrafficFlow");
    model.result("pg3").feature("mesh1").set("colortabletrans", "nonlinear");
    model.result("pg3").feature("mesh1").set("nonlinearcolortablerev", true);

//    In the Settings window for Mesh, click to expand the Title section.
//    From the Title type list, select None.

    model.result("pg3").feature("mesh1").set("titletype", "none");

//    Right-click Mesh 1 and choose Selection.

    model.result("pg3").feature("mesh1").create("sel1", "MeshSelection");

//    Select Boundaries 15, 24.

    model.result("pg3").feature("mesh1").feature("sel1").selection().set(15, 24);
    model.result("pg3").run();

//    In the Model Builder window, click Mesh 1.
//    In the Settings window for Mesh, locate the Coloring and Style section.
//    From the Element color list, select Gray.

    model.result("pg3").feature("mesh1").set("elemcolor", "gray");

//    In the Interface Health: Gauss Points toolbar, click Plot.

    model.result("pg3").run();

//    In the Home toolbar, click Add Plot Group and choose 1D Plot Group.

    model.result().create("pg4", "PlotGroup1D");
    model.result("pg4").run();

//    In the Settings window for 1D Plot Group, type Load Displacement Curve in the Label text field.

    model.result("pg4").label("Load Displacement Curve");

//    Click to expand the Title section.

    return model;
  }

  public static Model run2(Model model) {
//    From the Title type list, select Manual.

    model.result("pg4").set("titletype", "manual");

//    In the Title text area, type Load displacement curve.

    model.result("pg4").set("title", "Load displacement curve");

//    Right-click Load Displacement Curve and choose Global.

    model.result("pg4").create("glob1", "Global");
    model.result("pg4").feature("glob1").set("markerpos", "datapoints");
    model.result("pg4").feature("glob1").set("linewidth", "preference");

//    In the Settings window for Global, locate the y-Axis Data section.
//    In the table, enter the following settings:

    model.result("pg4").feature("glob1").setIndex("expr", "2*F_lp", 0);
    model.result("pg4").feature("glob1").setIndex("descr", "Load", 0);

//    The factor 2 is needed to compensate for the model symmetry.
//    Click to expand the Legends section.
//    Clear the Show legends checkbox.

    model.result("pg4").feature("glob1").set("legend", false);

//    Locate the x-Axis Data section.
//    From the Parameter list, select Expression.

    model.result("pg4").feature("glob1").set("xdata", "expr");

//    In the Expression text field, type u_lp.

    model.result("pg4").feature("glob1").set("xdataexpr", "u_lp");

//    In the Load Displacement Curve toolbar, click Plot.

    model.result("pg4").run();

//    Finally, evaluate the maximal load that this beam can carry under this loading condition.
//    In the Results toolbar, click Evaluation Group.

    model.result().evaluationGroup().create("eg1", "EvaluationGroup");

//    In the Settings window for Evaluation Group, type Load (solid) in the Label text field.

    model.result().evaluationGroup("eg1").label("Load (solid)");

//    In the Load (solid) toolbar, click Global Evaluation.

    model.result().evaluationGroup("eg1").create("gev1", "EvalGlobal");

//    In the Settings window for Global Evaluation, locate the Expressions section.
//    In the table, enter the following settings:

    model.result().evaluationGroup("eg1").feature("gev1").setIndex("expr", "2*F_lp", 0);
    model.result().evaluationGroup("eg1").feature("gev1").setIndex("descr", "Load", 0);

//    Again, the factor 2 is due to model symmetry.
//    Locate the Data Series Operation section.
//    From the Transformation list, select Maximum.

    model.result().evaluationGroup("eg1").feature("gev1").set("dataseries", "maximum");

//    In the Load (solid) toolbar, click Evaluate.

    model.result().evaluationGroup("eg1").run();
    model.result("pg1").run();

    model.title("Mixed-Mode Debonding of a Laminated Composite");

    model
         .description("Interfacial failure by delamination or debonding can be simulated with a Cohesive Zone Model (CZM). This example shows the implementation of a CZM with a bilinear traction-separation law. It is used to predict the mixed-mode softening onset and delamination propagation in a composite material.");

    return model;
  }

  public static void main(String[] args) {
    Model model = run();
    run2(model);
  }

}