/*
 * cycle_counting_benchmark.java
 */

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

/** Model exported on May 15 2026, 10:53 by COMSOL 6.4.0.421. */
public class cycle_counting_benchmark {

  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 2D.
//    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", 2);

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

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

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

//    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");

//    In the Geometry toolbar, click Rectangle.

    model.component("comp1").geom("geom1").create("r1", "Rectangle");

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

    model.component("comp1").geom("geom1").feature("r1").set("size", new double[]{10, 1});

//    In the Height text field, type 100.

    model.component("comp1").geom("geom1").feature("r1").set("size", new double[]{10, 100});

//    In the Geometry toolbar, click Rectangle.

    model.component("comp1").geom("geom1").run("r1");
    model.component("comp1").geom("geom1").create("r2", "Rectangle");

//    In the Settings window for Rectangle, locate the Size and Shape section.
//    In the Width text field, type 10-6.25.

    model.component("comp1").geom("geom1").feature("r2").set("size", new String[]{"10-6.25", "1"});

//    In the Height text field, type 50-sqrt(12.5^2-8.75^2).

    model.component("comp1").geom("geom1").feature("r2")
         .set("size", new String[]{"10-6.25", "50-sqrt(12.5^2-8.75^2)"});

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

    model.component("comp1").geom("geom1").feature("r2").set("pos", new double[]{6.25, 0});

//    In the Geometry toolbar, click Circle.

    model.component("comp1").geom("geom1").run("r2");
    model.component("comp1").geom("geom1").create("c1", "Circle");

//    In the Settings window for Circle, locate the Position section.
//    In the x text field, type 6.25+12.5.

    model.component("comp1").geom("geom1").feature("c1").set("pos", new String[]{"6.25+12.5", "0"});

//    In the y text field, type 50-sqrt(12.5^2-8.75^2).

    model.component("comp1").geom("geom1").feature("c1")
         .set("pos", new String[]{"6.25+12.5", "50-sqrt(12.5^2-8.75^2)"});

//    Locate the Size and Shape section.
//    In the Radius text field, type 12.5.

    model.component("comp1").geom("geom1").feature("c1").set("r", 12.5);

//    In the Geometry toolbar, click Booleans and Partitions and choose Difference.

    model.component("comp1").geom("geom1").run("c1");
    model.component("comp1").geom("geom1").create("dif1", "Difference");

//    From the bigger rectangle subtract the smaller rectangle and the circle.

    model.component("comp1").geom("geom1").feature("dif1").selection("input").set("r1");
    model.component("comp1").geom("geom1").feature("dif1").selection("input2").set("c1", "r2");
    model.component("comp1").geom("geom1").runPre("fin");

//    Specify the load cycle.
//    In the Home toolbar, click Functions and choose Global > Interpolation.

    model.func().create("int1", "Interpolation");

//    In the Settings window for Interpolation, locate the Definition section.
//    In the table, enter the following settings:

    model.func("int1").setIndex("table", 1, 0, 0);
    model.func("int1").setIndex("table", -2, 0, 1);
    model.func("int1").setIndex("table", 2, 1, 0);
    model.func("int1").setIndex("table", 1, 1, 1);
    model.func("int1").setIndex("table", 3, 2, 0);
    model.func("int1").setIndex("table", -3, 2, 1);
    model.func("int1").setIndex("table", 4, 3, 0);
    model.func("int1").setIndex("table", 5, 3, 1);
    model.func("int1").setIndex("table", 5, 4, 0);
    model.func("int1").setIndex("table", -1, 4, 1);
    model.func("int1").setIndex("table", 6, 5, 0);
    model.func("int1").setIndex("table", 3, 5, 1);
    model.func("int1").setIndex("table", 7, 6, 0);
    model.func("int1").setIndex("table", -4, 6, 1);
    model.func("int1").setIndex("table", 8, 7, 0);
    model.func("int1").setIndex("table", 4, 7, 1);
    model.func("int1").setIndex("table", 9, 8, 0);
    model.func("int1").setIndex("table", -2, 8, 1);

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

//    In the Model Builder window, under Component 1 (comp1), click Solid Mechanics (solid).
//    In the Settings window for Solid Mechanics, locate the 2D Approximation section.
//    From the list, select Plane stress.

    model.component("comp1").physics("solid").prop("Type2D").set("Type2D", "PlaneStress");

//    Locate the Thickness section.
//    In the \[d\] text field, type 0.00625.

    model.component("comp1").physics("solid").prop("d").set("d", 0.00625);

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

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

//    Select Boundaries 1, 2.

    model.component("comp1").physics("solid").feature("sym1").selection().set(1, 2);

//    In the Physics toolbar, click Boundaries and choose Boundary Load.

    model.component("comp1").physics("solid").create("bndl1", "BoundaryLoad", 1);

//    Select Boundary 3.

    model.component("comp1").physics("solid").feature("bndl1").selection().set(3);

//    In the Settings window for Boundary Load, locate the Force section.
//    From the Load type list, select Total force.

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

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

    model.component("comp1").physics("solid").feature("bndl1").set("force", new String[]{"0", "F*int1(case)", "0"});

//    In the Model Builder window, under Global Definitions, click Parameters 1.
//    In the Settings window for Parameters, locate the Parameters section.
//    In the table, enter the following settings:

    model.param().set("F", "10*6.25*12.5/2");
    model.param().descr("F", "Load unit");
    model.param().set("case", "1");
    model.param().descr("case", "Load case");

//    In the Model Builder window, right-click Component 1 (comp1) > Materials and choose Blank Material.

    model.component("comp1").material().create("mat1", "Common");

//    In the Settings window for Material, locate the Material Contents section.
//    In the table, enter the following settings:

    model.component("comp1").material("mat1").propertyGroup()
         .create("Enu", "Enu", "Young's_modulus_and_Poisson's_ratio");
    model.component("comp1").material("mat1").propertyGroup("Enu").set("E", new String[]{"69e9"});
    model.component("comp1").material("mat1").propertyGroup("Enu").set("nu", new String[]{"0.34"});
    model.component("comp1").material("mat1").propertyGroup("def").set("density", new String[]{"2700"});

//    In the Model Builder window, under Study 1, click Step 1: Stationary.
//    In the Settings window for Stationary, 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", "F", 0);
    model.study("std1").feature("stat").setIndex("plistarr", "", 0);
    model.study("std1").feature("stat").setIndex("punit", "", 0);
    model.study("std1").feature("stat").setIndex("pname", "F", 0);
    model.study("std1").feature("stat").setIndex("plistarr", "", 0);
    model.study("std1").feature("stat").setIndex("punit", "", 0);

//    In the table, enter the following settings:

    model.study("std1").feature("stat").setIndex("pname", "case", 0);
    model.study("std1").feature("stat").setIndex("plistarr", "range(1,1,9)", 0);

//    In the Study toolbar, click Compute.

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

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

    model.result().create("pg1", "PlotGroup2D");
    model.result("pg1").set("data", "dset1");
    model.result("pg1").setIndex("looplevel", 9, 0);
    model.result("pg1").label("Stress (solid)");
    model.result("pg1").set("frametype", "spatial");
    model.result("pg1").create("surf1", "Surface");
    model.result("pg1").feature("surf1").set("expr", new String[]{"solid.misesGp"});
    model.result("pg1").feature("surf1").set("threshold", "manual");
    model.result("pg1").feature("surf1").set("thresholdvalue", 0.2);
    model.result("pg1").feature("surf1").set("colortable", "Rainbow");
    model.result("pg1").feature("surf1").set("colortabletrans", "none");
    model.result("pg1").feature("surf1").set("colorscalemode", "linear");
    model.result("pg1").feature("surf1").set("resolution", "normal");
    model.result("pg1").feature("surf1").set("refine", 2);
    model.result("pg1").feature("surf1").set("colortable", "Prism");
    model.result("pg1").feature("surf1").create("def", "Deform");
    model.result("pg1").feature("surf1").feature("def").set("expr", new String[]{"u", "v"});
    model.result("pg1").feature("surf1").feature("def").set("descr", "Displacement field");
    model.result("pg1").run();

//    Set default units for result presentation.
//    In the Results toolbar, click Configurations and choose Preferred Units.

    model.result().configuration().create("prfu1", "PreferredUnits");

//    In the Settings window for Preferred Units, locate the Units section.
//    Click Add Physical Quantity.
//    In the Physical Quantity dialog, select Solid Mechanics > Stress tensor (N/m^2) in the tree.
//    Click OK.

    model.result().configuration("prfu1")
         .setIndex("quantityunits", new String[]{"stress", "Stress tensor", "N/m^2", "N/m^2"}, 0);

//    In the Settings window for Preferred Units, locate the Units section.
//    In the table, enter the following settings:

    model.result().configuration("prfu1").setIndex("quantityunits", "MPa", 0, 3);

//    Click Apply.

    model.result().configuration("prfu1").apply();
    model.result("pg1").run();

//    Mirror solution of a quarter of a specimen and create results for a full specimen.
//    In the Results toolbar, click More Datasets and choose Mirror 2D.

    model.result().dataset().create("mir1", "Mirror2D");

//    In the Settings window for Mirror 2D, locate the Axis Data section.
//    In row Point 2, set Y to 100.

    model.result().dataset("mir1").setIndex("genpoints", 100, 1, 1);

//    In the Results toolbar, click More Datasets and choose Mirror 2D.

    model.result().dataset().create("mir2", "Mirror2D");

//    In the Settings window for Mirror 2D, locate the Data section.
//    From the Dataset list, select Mirror 2D 1.

    model.result().dataset("mir2").set("data", "mir1");

//    Locate the Axis Data section.
//    In row Point 1, set x to -6.25.

    model.result().dataset("mir2").setIndex("genpoints", -6.25, 0, 0);

//    In row Point 2, set x to 6.25.

    model.result().dataset("mir2").setIndex("genpoints", 6.25, 1, 0);

//     and <l>y</l> to <c>0</c>.

    model.result().dataset("mir2").setIndex("genpoints", 0, 1, 1);

//    Display stress state in the whole specimen.

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

//    In the Model Builder window, under Results, click Stress (solid).
//    In the Settings window for 2D Plot Group, locate the Data section.
//    From the Dataset list, select Mirror 2D 2.

    model.result("pg1").set("data", "mir2");

//    From the Parameter value (case) list, select 2.

    model.result("pg1").setIndex("looplevel", 2, 0);

//    Click to expand the Title section.
//    From the Title type list, select None.

    model.result("pg1").set("titletype", "none");

//    Locate the Color Legend section.
//    Select the Show units checkbox.

    model.result("pg1").set("showlegendsunit", true);
    model.result("pg1").run();

//    In the Model Builder window, expand the Stress (solid) node, then click Surface 1.
//    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 > Stress > Second Piola–Kirchhoff stress (material and geometry frames) - N/m² > solid.SGpYY - Second Piola–Kirchhoff stress, YY-component.

    model.result("pg1").feature("surf1").set("expr", "solid.SGpYY");
    model.result("pg1").feature("surf1").set("descr", "Second Piola\u2013Kirchhoff stress, YY-component");

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

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

//    In the Results toolbar, click 1D Plot Group.

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

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

    model.result("pg2").label("Load Cycle Response");

//    Right-click Load Cycle Response and choose Point Graph.

    model.result("pg2").create("ptgr1", "PointGraph");
    model.result("pg2").feature("ptgr1").set("markerpos", "datapoints");
    model.result("pg2").feature("ptgr1").set("linewidth", "preference");

//    Select Point 3.

    model.result("pg2").feature("ptgr1").selection().set(3);

//    In the Settings window for Point Graph, click Replace Expression in the upper-right corner of the y-Axis Data section.
//    From the menu, choose Component 1 (comp1) > Solid Mechanics > Stress > Second Piola–Kirchhoff stress (material and geometry frames) - N/m² > solid.SGpYY - Second Piola–Kirchhoff stress, YY-component.

    model.result("pg2").feature("ptgr1").set("expr", "solid.SGpYY");
    model.result("pg2").feature("ptgr1").set("descr", "Second Piola\u2013Kirchhoff stress, YY-component");

//    In the Load Cycle Response toolbar, click Plot.

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

//    Specify the load cycle.
//    In the Home toolbar, click Functions and choose Global > Analytic.

    model.func().create("an1", "Analytic");

//    In the Settings window for Analytic, locate the Definition section.
//    In the Arguments text field, type R, N.

    model.func("an1").set("args", "R, N");

//    In the Expression text field, type (94e6*(R/-0.36)^1.15)*N^-0.119.

    model.func("an1").set("expr", "(94e6*(R/-0.36)^1.15)*N^-0.119");

//    In the Home toolbar, click Add Physics to open the Add Physics window.
//    In the tree, select Structural Mechanics > Fatigue (ftg).
//    Find the Physics interfaces in study subsection.
//    In the table, enter the following settings:
//    Click Add to Component 1 in the window toolbar.

    model.component("comp1").physics().create("ftg", "Fatigue", "geom1");

    model.study("std1").feature("stat").setSolveFor("/physics/ftg", false);

//    In the Home toolbar, click Add Physics to close the Add Physics window.
//    In the Physics toolbar, click Points and choose Cumulative Damage.

    model.component("comp1").physics("ftg").create("cdam1", "CumulativeDamageModel2", 0);

//    Any point in the thin section away from the notch will give the same fatigue response.
//    Select Point 3.

    model.component("comp1").physics("ftg").feature("cdam1").selection().set(3);

//    In the Settings window for Cumulative Damage, locate the Solution Field section.
//    From the Physics interface list, select Solid Mechanics (solid).

    model.component("comp1").physics("ftg").feature("cdam1").set("fatigueInputPhysics", "solid");

//    Locate the Cycle Counting Parameters section.
//    Find the Discretization subsection.
//    In the \[N_{\textrm{m}}\] text field, type 5.

    model.component("comp1").physics("ftg").feature("cdam1").set("cDamDiscretizationMean", 5);

//    In the \[N_{\textrm{r}}\] text field, type 7.

    model.component("comp1").physics("ftg").feature("cdam1").set("cDamDiscretizationRange", 7);

//    Locate the Damage Model Parameters section.
//    From the \[f_\textrm{SN}(R,\,N)\] list, select an1.

    model.component("comp1").physics("ftg").feature("cdam1").set("lifeCurve", "an1");

//    In the \[m\] text field, type 100000.

    model.component("comp1").physics("ftg").feature("cdam1").set("m", 100000);

//    Set cutoff which can be seen as a limit for infinite life that does not contribute to the damage.
//    Find the Evaluation settings subsection.
//    In the \[N_{\textrm{cut}}\] text field, type 1e8.

    model.component("comp1").physics("ftg").feature("cdam1").set("Ncut", "1e8");

//    In the Home toolbar, click Add Study to open the Add Study window.
//    Find the Physics interfaces in study subsection.
//    In the table, enter the following settings:
//    Find the Studies subsection.
//    In the Select Study tree, select Preset Studies for Selected Physics Interfaces > Fatigue.
//    Click Add Study in the window toolbar.

    model.study().create("std2");
    model.study("std2").create("ftge", "Fatigue");
    model.study("std2").feature("ftge").set("solnum", "auto");
    model.study("std2").feature("ftge").set("usesol", "off");
    model.study("std2").feature("ftge").set("outputmap", new String[]{});
    model.study("std2").feature("ftge").setSolveFor("/physics/solid", false);

//    In the Home toolbar, click Add Study to close the Add Study window.
//    In the Settings window for Fatigue, locate the Values of Dependent Variables section.
//    Find the Values of variables not solved for subsection.
//    From the Settings list, select User controlled.

    model.study("std2").feature("ftge").set("usesol", true);

//    From the Method list, select Solution.

    model.study("std2").feature("ftge").set("notsolmethod", "sol");

//    From the Study list, select Study 1, Stationary.

    model.study("std2").feature("ftge").set("notstudy", "std1");

//    In the Study toolbar, click Compute.

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

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

    model.result().create("pg3", "PlotGroup1D");
    model.result("pg3").set("data", "dset2");
    model.result("pg3").create("ptgr1", "PointGraph");
    model.result("pg3").feature("ptgr1").selection().all();
    model.result("pg3").feature("ptgr1").set("expr", new String[]{"ftg.fus"});
    model.result("pg3").feature("ptgr1").selection().set(3);
    model.result("pg3").label("Fatigue Usage Factor (ftg)");
    model.result().create("pg4", "PlotGroup2D");
    model.result("pg4").set("data", "none");
    model.result("pg4").create("hist1", "MatrixHistogram");
    model.result("pg4").feature("hist1").set("data", "dset2");
    model.result("pg4").feature("hist1").set("expr", "ftg.cdam1.csc");
    model.result("pg4").label("Stress Cycle Distribution (ftg)");
    model.result("pg4").feature("hist1").set("xdata", "Mean stress");
    model.result("pg4").feature("hist1").set("ydata", "Stress amplitude");
    model.result("pg4").feature("hist1").create("hght", "HistogramHeight");
    model.result("pg4").feature("hist1").feature("hght").label("Height expression (for 3D histogram)");
    model.result().create("pg5", "PlotGroup2D");
    model.result("pg5").set("data", "none");
    model.result("pg5").create("hist1", "MatrixHistogram");
    model.result("pg5").feature("hist1").set("data", "dset2");
    model.result("pg5").feature("hist1").set("expr", "ftg.cdam1.rus");
    model.result("pg5").label("Fatigue Usage Distribution (ftg)");
    model.result("pg5").feature("hist1").set("xdata", "Mean stress");
    model.result("pg5").feature("hist1").set("ydata", "Stress amplitude");
    model.result("pg5").feature("hist1").create("hght", "HistogramHeight");
    model.result("pg5").feature("hist1").feature("hght").label("Height expression (for 3D histogram)");
    model.result("pg3").run();
    model.result("pg4").run();
    model.result("pg4").run();

//    In the Model Builder window, expand the Stress Cycle Distribution (ftg) node, then click Matrix Histogram 1.
//    In the Settings window for Matrix Histogram, locate the Axes section.
//    From the Unit list, select MPa.

    model.result("pg4").feature("hist1").set("axisunit", "MPa");

//    In the Stress Cycle Distribution (ftg) toolbar, click Plot.

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

//    Click the Zoom Extents button in the Graphics toolbar.

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

//    In the Model Builder window, expand the Fatigue Usage Distribution (ftg) node, then click Matrix Histogram 1.
//    In the Settings window for Matrix Histogram, locate the Axes section.
//    From the Unit list, select MPa.

    model.result("pg5").feature("hist1").set("axisunit", "MPa");

//    In the Fatigue Usage Distribution (ftg) toolbar, click Plot.

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

//    Click the Zoom Extents button in the Graphics toolbar.

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

    model.title("Cycle Counting in Fatigue Analysis \u2014 Benchmark");

    model
         .description("A benchmark problem of the Rainflow Counting algorithm compares results between ASTM and the COMSOL Fatigue Module using a flat tensile test specimen. An extension is made for the cumulative damage calculation following the Palmgren\u2013Miner model and results are compared with analytical expressions.");

    return model;
  }

  public static void main(String[] args) {
    run();
  }

}