diff --git a/docs/source/examples/Example-Beam_Optimization.rst b/docs/source/examples/Example-Beam_Optimization.rst index c24b0d2c1..c9b825b70 100644 --- a/docs/source/examples/Example-Beam_Optimization.rst +++ b/docs/source/examples/Example-Beam_Optimization.rst @@ -76,7 +76,7 @@ assign a design variable number for the thickness parameter, and return a :class # Set one thickness dv for every property group con = constitutive.IsoRectangleBeamConstitutive(prop, t=t, w=w, tNum=dvNum) - # Defines local y/width direction for beam + # Defines local y/thickness direction for beam refAxis = np.array([0.0, 1.0, 0.0]) transform = elements.BeamRefAxisTransform(refAxis) @@ -203,12 +203,16 @@ Finally, we can plot the optimized thickness distribution using matplotlib and c # Get analytical solution t_exact = np.sqrt(6 * (L - x) * V / w / ys) + # Compute max thickness value + t0 = np.sqrt(6 * L * V / w / ys) + # Plot results for solution - plt.plot(x, t_opt, "o", x, t_exact) + plt.plot(x / L, t_opt / t0, "o", x, t_exact / t0) plt.legend(["optimized", "analytical"]) - plt.ylabel("t(x)") - plt.xlabel("x") + plt.ylabel(r"$\frac{t(x)}{t_0}$", fontsize=16) + plt.xlabel(r"$\frac{x}{L}$", fontsize=16, labelpad=-5) plt.title("Optimal beam thickness profile") + plt.text(0.05, 0.25, r"$t_0 = \sqrt{\frac{6VL}{w\sigma_y}}$", fontsize=12) plt.show() .. image:: images/beam_plot.png diff --git a/docs/source/examples/Example-Composite_Optimization.rst b/docs/source/examples/Example-Composite_Optimization.rst new file mode 100644 index 000000000..eef24b6d1 --- /dev/null +++ b/docs/source/examples/Example-Composite_Optimization.rst @@ -0,0 +1,245 @@ +Composite plate optimization with MPhys +*************************************** +.. note:: The script for this example can be found under the `examples/plate/` directory. + +This example further demonstrates TACS structural optimization capabilities through :ref:`mphys/mphys:MPhys`. +The beam model that we will be using for this problem is a rectangular beam, +cantilevered, with a shear load applied at the tip. The beam is discretized using +100 beam elements along it's span. + +In this example a composite plate compliance minimization problem, based off of one originally proposed by Lund and Stegmann [`1`_], is solved. +A 1 m x 1 m x 0.05 m composite plate is clamped on all edges and subjected to +a uniform pressure of 100 kPa loading on the top. +The plate is discretized into 100 quad shell elements each having its own independent laminate layup. +The goal of this example will be to optimize the layup of each element in the plate in order to minimize the total compliance energy. +A diagram of the problem is shown below. + +.. image:: images/plate_pressure.png + :width: 800 + :alt: Plate problem + +To make this problem tractable for gradient-based optimization we will make the following simplifications: + + 1. For each layup we can ony select plies of the following angles: :math:`0^\circ, 45^\circ, -45^\circ, 90^\circ`. + + 2. We will neglect the discrete nature of the plies the effects of stacking sequence and instead homogenize or "smear" each ply angles contribution to laminate stiffness based on its proportion of plies in the layup or "ply fraction". + +The optimization problem can now be summarized as follows: + +.. image:: images/comp_opt_def.png + :width: 800 + :alt: Optimization problem + +First, we import required libraries, define the model bdf file, and define important problem constants: + +.. code-block:: python + + import os + + import openmdao.api as om + import numpy as np + from mphys import Multipoint + from mphys.scenario_structural import ScenarioStructural + + from tacs import elements, constitutive, functions + from tacs.mphys import TacsBuilder + + # BDF file containing mesh + bdf_file = os.path.join(os.path.dirname(__file__), "partitioned_plate.bdf") + + # Material properties + rho = 1550.0 + E1 = 54e9 + E2 = 18e9 + nu12 = 0.25 + G12 = 9e9 + G13 = 9e9 + Xt = 2410.0e6 + Xc = 1040.0e6 + Yt = 73.0e6 + Yc = 173.0e6 + S12 = 71.0e6 + + # Shell thickness + ply_thickness = 1.25e-3 # m + plate_thickness = 0.05 # m + tMin = 0.002 # m + tMax = 0.05 # m + + # Ply angles/initial ply fractions + ply_angles = np.deg2rad([0.0, 45.0, -45.0, 90.0]) + ply_fractions = np.array([0.25, 0.25, 0.25, 0.25]) + + # Pressure load to apply to plate + P = 100e3 + +Next we define an :func:`~tacs.pytacs.elemCallBack` function for setting up the TACS elements and design variables. +We use the :class:`~tacs.constitutive.SmearedCompositeShellConstitutive` class here for the constitutive properties, and +assign four design variable numbers to each element (one for each ply fraction), and return a :class:`~tacs.elements.Quad4Shell` element class. + +.. code-block:: python + + # Callback function used to setup TACS element objects and DVs + def element_callback(dvNum, compID, compDescript, elemDescripts, specialDVs, **kwargs): + # Create ply object + ortho_prop = constitutive.MaterialProperties( + rho=rho, + E1=E1, + E2=E2, + nu12=nu12, + G12=G12, + G13=G13, + G23=G13, + Xt=Xt, + Xc=Xc, + Yt=Yt, + Yc=Yc, + S12=S12, + ) + ortho_ply = constitutive.OrthotropicPly(ply_thickness, ortho_prop) + # Create the layup list (one for each angle) + ortho_layup = [ortho_ply, ortho_ply, ortho_ply, ortho_ply] + # Assign each ply fraction a unique DV + ply_fraction_dv_nums = np.array( + [dvNum, dvNum + 1, dvNum + 2, dvNum + 3], dtype=np.intc + ) + # Create smeared stiffness object based on ply angles/fractions + con = constitutive.SmearedCompositeShellConstitutive( + ortho_layup, + plate_thickness, + ply_angles, + ply_fractions, + ply_fraction_dv_nums=ply_fraction_dv_nums, + ) + + # Define reference axis to define local 0 deg direction + refAxis = np.array([1.0, 0.0, 0.0]) + transform = elements.ShellRefAxisTransform(refAxis) + + # Pass back the appropriate tacs element object + elem = elements.Quad4Shell(transform, con) + + return elem + +We define a :func:`problem_setup` to add fixed loads and eval functions. +Here we specify the plate compliance energy (:class:`~tacs.functions.Compliance`) as an output for our analysis +and add our 100 kPa pressure load. + +.. code-block:: python + + def problem_setup(scenario_name, fea_assembler, problem): + """ + Helper function to add fixed forces and eval functions + to structural problems used in tacs builder + """ + + # Add TACS Functions + problem.addFunction("compliance", functions.Compliance) + + # Add forces to static problem + allComponents = fea_assembler.selectCompIDs() + problem.addPressureToComponents(allComponents, P) + +For our last helper function we define a :func:`constraint_setup` function. +This function can be used to add additional relational constraints to the design variables we defined in the ``element_callback``. +In particular, we want to enforce a new constraint (100 in total) such that the ply fractions within each element should sum to unity. +We can accomplish this by utilizing the :class:`~tacs.constraints.DVConstraint`. + +.. code-block:: python + + def constraint_setup(scenario_name, fea_assembler, constraint_list): + """ + Helper function to setup tacs constraint classes + """ + constr = fea_assembler.createDVConstraint("ply_fractions") + allComponents = fea_assembler.selectCompIDs() + constr.addConstraint( + "sum", allComponents, dvIndices=[0, 1, 2, 3], dvWeights=[1.0, 1.0, 1.0, 1.0] + ) + constraint_list.append(constr) + +Here we define our :class:`~mphys.Multipoint` (essentially an OpenMDAO ``Group``) which will contain our analysis :class:`~mphys.Scenario`. +To do this, we instantiate the :class:`~tacs.mphys.builder.TacsBuilder` using the ``element_callback`` and ``problem_setup`` we defined above. +We create OpenMDAO ``Component``'s to feed design variable and mesh inputs to the ``Scenario`` component. +We use this builder to create an MPhys :class:`~mphys.StructuralScenario`. + +.. code-block:: python + + class PlateModel(Multipoint): + def setup(self): + struct_builder = TacsBuilder( + mesh_file=bdf_file, + element_callback=element_callback, + problem_setup=problem_setup, + constraint_setup=constraint_setup, + coupled=False, + check_partials=True, + ) + struct_builder.initialize(self.comm) + dv_array = struct_builder.get_initial_dvs() + + dvs = self.add_subsystem("dvs", om.IndepVarComp(), promotes=["*"]) + dvs.add_output("dv_struct", dv_array) + + self.add_subsystem("mesh", struct_builder.get_mesh_coordinate_subsystem()) + self.mphys_add_scenario( + "pressure_load", ScenarioStructural(struct_builder=struct_builder) + ) + self.mphys_connect_scenario_coordinate_source("mesh", "pressure_load", "struct") + + self.connect("dv_struct", "pressure_load.dv_struct") + +At this point we setup the OpenMDAO ``Problem`` class that we will use to perform our optimization. +We assign our ``PlateModel`` to the problem class and set ``ScipyOptimizeDriver``. +We define our design variables, constraint, and objective. +Finally we run the problem driver to optimize the problem. + +.. code-block:: python + + prob = om.Problem() + prob.model = PlateModel() + model = prob.model + + # Declare design variables, objective, and constraint + model.add_design_var("dv_struct", lower=0.0, upper=1.0) + model.add_objective("pressure_load.compliance", scaler=1.0) + model.add_constraint("pressure_load.ply_fractions.sum", equals=1.0, linear=True) + + # Configure optimizer + prob.driver = om.ScipyOptimizeDriver(debug_print=["objs", "nl_cons"], maxiter=100) + prob.driver.options["optimizer"] = "SLSQP" + + # Setup OpenMDAO problem + prob.setup() + + # Output N2 representation of OpenMDAO model + om.n2(prob, show_browser=False, outfile="tacs_struct.html") + + # Run optimization + prob.run_driver() + +After the optimization completes the user should see a print out to screen like shown below. + +>>> Optimization terminated successfully (Exit mode 0) +>>> Current function value: 8.571649588963465 +>>> Iterations: 34 +>>> Function evaluations: 34 +>>> Gradient evaluations: 34 +>>> Optimization Complete +>>> ----------------------------------- + +Once the optimization is complete we can post-process results. +We can write our optimized beam model to a BDF file so they can +be processed in other commonly used FEM software. +The ``f5`` solution file at each optimization iteration can also be converted to a Tecplot or Paraview files using ``f5totec`` or ``f5tovtk``, respectively. +The optimized ply fraction distributions for each angle can be visualized by plotting the contours of the following variables in Tecplot or Paraview: ``dv2``, ``dv3``, ``dv4``, ``dv5``. +A visualization of the optimized result is shown below: + +.. image:: images/pf_opt.png + :width: 800 + :alt: Plate solution + +.. rubric:: References + +.. [1] Lund, E. and Stegmann, J., β€œOn structural optimization of composite shell structures using a discrete constitutive parametrization,” Wind Energy, Vol. 8, No. 1, 2005, pp. 109–124. + diff --git a/docs/source/examples/images/comp_opt_def.png b/docs/source/examples/images/comp_opt_def.png new file mode 100644 index 000000000..ce092494b Binary files /dev/null and b/docs/source/examples/images/comp_opt_def.png differ diff --git a/docs/source/examples/images/pf_opt.png b/docs/source/examples/images/pf_opt.png new file mode 100644 index 000000000..feef272ea Binary files /dev/null and b/docs/source/examples/images/pf_opt.png differ diff --git a/docs/source/examples/images/plate_pressure.png b/docs/source/examples/images/plate_pressure.png new file mode 100644 index 000000000..eb32f1871 Binary files /dev/null and b/docs/source/examples/images/plate_pressure.png differ diff --git a/docs/source/index.rst b/docs/source/index.rst index 1c342e5d0..2c76a2b21 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -26,6 +26,7 @@ Examples examples/Example-Plate examples/Example-Transient_Battery examples/Example-Beam_Optimization + examples/Example-Composite_Optimization References ========== diff --git a/docs/source/mphys/builder.rst b/docs/source/mphys/builder.rst index 4f24789f6..f6ce38454 100644 --- a/docs/source/mphys/builder.rst +++ b/docs/source/mphys/builder.rst @@ -3,7 +3,6 @@ TacsBuilder class TACS MPhys interface centers around the :class:`~tacs.mphys.builder.TacsBuilder` class. The :class:`~tacs.mphys.builder.TacsBuilder` is responsible for creating :ref:`Scenarios `, an OpenMDAO group that contains an analysis condition in a multipoint optimization. -:ref:`mphys:builders` . API Reference ------------- diff --git a/examples/beam/beam_opt.py b/examples/beam/beam_opt.py index 39a721547..e7fbe3106 100644 --- a/examples/beam/beam_opt.py +++ b/examples/beam/beam_opt.py @@ -56,7 +56,7 @@ def element_callback(dvNum, compID, compDescript, elemDescripts, specialDVs, **k # Set one thickness dv for every property group con = constitutive.IsoRectangleBeamConstitutive(prop, t=t, w=w, tNum=dvNum) - # Defines local y/width direction for beam + # Defines local y/thickness direction for beam refAxis = np.array([0.0, 1.0, 0.0]) transform = elements.BeamRefAxisTransform(refAxis) diff --git a/examples/plate/optimize_composite_plate.py b/examples/plate/optimize_composite_plate.py new file mode 100644 index 000000000..4f00ff69f --- /dev/null +++ b/examples/plate/optimize_composite_plate.py @@ -0,0 +1,162 @@ +""" +Mass minimization of uCRM wingbox subject to a constant vertical force +""" +import os + +import openmdao.api as om +import numpy as np +from mphys import Multipoint +from mphys.scenario_structural import ScenarioStructural + +from tacs import elements, constitutive, functions +from tacs.mphys import TacsBuilder + +# BDF file containing mesh +bdf_file = os.path.join(os.path.dirname(__file__), "partitioned_plate.bdf") + +# Material properties +rho = 1550.0 +E1 = 54e9 +E2 = 18e9 +nu12 = 0.25 +G12 = 9e9 +G13 = 9e9 +Xt = 2410.0e6 +Xc = 1040.0e6 +Yt = 73.0e6 +Yc = 173.0e6 +S12 = 71.0e6 + +# Shell thickness +ply_thickness = 1.25e-3 # m +plate_thickness = 0.05 # m +tMin = 0.002 # m +tMax = 0.05 # m + +# Ply angles/initial ply fractions +ply_angles = np.deg2rad([0.0, 45.0, -45.0, 90.0]) +ply_fractions = np.array([0.25, 0.25, 0.25, 0.25]) + +# Pressure load to apply to plate +P = 100e3 + + +# Callback function used to setup TACS element objects and DVs +def element_callback(dvNum, compID, compDescript, elemDescripts, specialDVs, **kwargs): + # Create ply object + ortho_prop = constitutive.MaterialProperties( + rho=rho, + E1=E1, + E2=E2, + nu12=nu12, + G12=G12, + G13=G13, + G23=G13, + Xt=Xt, + Xc=Xc, + Yt=Yt, + Yc=Yc, + S12=S12, + ) + ortho_ply = constitutive.OrthotropicPly(ply_thickness, ortho_prop) + # Create the layup list (one for each angle) + ortho_layup = [ortho_ply, ortho_ply, ortho_ply, ortho_ply] + # Assign each ply fraction a unique DV + ply_fraction_dv_nums = np.array( + [dvNum, dvNum + 1, dvNum + 2, dvNum + 3], dtype=np.intc + ) + # Create smeared stiffness object based on ply angles/fractions + con = constitutive.SmearedCompositeShellConstitutive( + ortho_layup, + plate_thickness, + ply_angles, + ply_fractions, + ply_fraction_dv_nums=ply_fraction_dv_nums, + ) + + # Define reference axis to define local 0 deg direction + refAxis = np.array([1.0, 0.0, 0.0]) + transform = elements.ShellRefAxisTransform(refAxis) + + # Pass back the appropriate tacs element object + elem = elements.Quad4Shell(transform, con) + + return elem + + +def problem_setup(scenario_name, fea_assembler, problem): + """ + Helper function to add fixed forces and eval functions + to structural problems used in tacs builder + """ + + # Add TACS Functions + problem.addFunction("compliance", functions.Compliance) + + # Add forces to static problem + allComponents = fea_assembler.selectCompIDs() + problem.addPressureToComponents(allComponents, P) + + +def constraint_setup(scenario_name, fea_assembler, constraint_list): + """ + Helper function to setup tacs constraint classes + """ + constr = fea_assembler.createDVConstraint("ply_fractions") + allComponents = fea_assembler.selectCompIDs() + constr.addConstraint( + "sum", allComponents, dvIndices=[0, 1, 2, 3], dvWeights=[1.0, 1.0, 1.0, 1.0] + ) + constraint_list.append(constr) + + +class PlateModel(Multipoint): + def setup(self): + struct_builder = TacsBuilder( + mesh_file=bdf_file, + element_callback=element_callback, + problem_setup=problem_setup, + constraint_setup=constraint_setup, + coupled=False, + check_partials=True, + ) + struct_builder.initialize(self.comm) + dv_array = struct_builder.get_initial_dvs() + + dvs = self.add_subsystem("dvs", om.IndepVarComp(), promotes=["*"]) + dvs.add_output("dv_struct", dv_array) + + self.add_subsystem("mesh", struct_builder.get_mesh_coordinate_subsystem()) + self.mphys_add_scenario( + "pressure_load", ScenarioStructural(struct_builder=struct_builder) + ) + self.mphys_connect_scenario_coordinate_source("mesh", "pressure_load", "struct") + + self.connect("dv_struct", "pressure_load.dv_struct") + + +################################################################################ +# OpenMDAO setup +################################################################################ + +prob = om.Problem() +prob.model = PlateModel() +model = prob.model + +# Declare design variables, objective, and constraint +model.add_design_var("dv_struct", lower=0.0, upper=1.0) +model.add_objective("pressure_load.compliance", scaler=1.0) +model.add_constraint("pressure_load.ply_fractions.sum", equals=1.0, linear=True) + +# Configure optimizer +prob.driver = om.ScipyOptimizeDriver(debug_print=["objs", "nl_cons"], maxiter=100) +prob.driver.options["optimizer"] = "SLSQP" + +# Setup OpenMDAO problem +prob.setup() + +# Output N2 representation of OpenMDAO model +om.n2(prob, show_browser=False, outfile="tacs_struct.html") + +# Run optimization +prob.run_driver() diff --git a/examples/plate/partitioned_plate.bdf b/examples/plate/partitioned_plate.bdf new file mode 100644 index 000000000..2b07f5143 --- /dev/null +++ b/examples/plate/partitioned_plate.bdf @@ -0,0 +1,251 @@ +$pyNastran: version=msc +$pyNastran: punch=True +$pyNastran: encoding=utf-8 +$pyNastran: nnodes=121 +$pyNastran: nelements=100 +$NODES +GRID 1 0. 0. 0. +GRID 2 0. .1 0. +GRID 3 0. .2 0. +GRID 4 0. .3 0. +GRID 5 0. .4 0. +GRID 6 0. .5 0. +GRID 7 0. .6 0. +GRID 8 0. .7 0. +GRID 9 0. .8 0. +GRID 10 0. .9 0. +GRID 11 0. 1. 0. +GRID 12 .1 0. 0. +GRID 13 .1 .1 0. +GRID 14 .1 .2 0. +GRID 15 .1 .3 0. +GRID 16 .1 .4 0. +GRID 17 .1 .5 0. +GRID 18 .1 .6 0. +GRID 19 .1 .7 0. +GRID 20 .1 .8 0. +GRID 21 .1 .9 0. +GRID 22 .1 1. 0. +GRID 23 .2 0. 0. +GRID 24 .2 .1 0. +GRID 25 .2 .2 0. +GRID 26 .2 .3 0. +GRID 27 .2 .4 0. +GRID 28 .2 .5 0. +GRID 29 .2 .6 0. +GRID 30 .2 .7 0. +GRID 31 .2 .8 0. +GRID 32 .2 .9 0. +GRID 33 .2 1. 0. +GRID 34 .3 0. 0. +GRID 35 .3 .1 0. +GRID 36 .3 .2 0. +GRID 37 .3 .3 0. +GRID 38 .3 .4 0. +GRID 39 .3 .5 0. +GRID 40 .3 .6 0. +GRID 41 .3 .7 0. +GRID 42 .3 .8 0. +GRID 43 .3 .9 0. +GRID 44 .3 1. 0. +GRID 45 .4 0. 0. +GRID 46 .4 .1 0. +GRID 47 .4 .2 0. +GRID 48 .4 .3 0. +GRID 49 .4 .4 0. +GRID 50 .4 .5 0. +GRID 51 .4 .6 0. +GRID 52 .4 .7 0. +GRID 53 .4 .8 0. +GRID 54 .4 .9 0. +GRID 55 .4 1. 0. +GRID 56 .5 0. 0. +GRID 57 .5 .1 0. +GRID 58 .5 .2 0. +GRID 59 .5 .3 0. +GRID 60 .5 .4 0. +GRID 61 .5 .5 0. +GRID 62 .5 .6 0. +GRID 63 .5 .7 0. +GRID 64 .5 .8 0. +GRID 65 .5 .9 0. +GRID 66 .5 1. 0. +GRID 67 .6 0. 0. +GRID 68 .6 .1 0. +GRID 69 .6 .2 0. +GRID 70 .6 .3 0. +GRID 71 .6 .4 0. +GRID 72 .6 .5 0. +GRID 73 .6 .6 0. +GRID 74 .6 .7 0. +GRID 75 .6 .8 0. +GRID 76 .6 .9 0. +GRID 77 .6 1. 0. +GRID 78 .7 0. 0. +GRID 79 .7 .1 0. +GRID 80 .7 .2 0. +GRID 81 .7 .3 0. +GRID 82 .7 .4 0. +GRID 83 .7 .5 0. +GRID 84 .7 .6 0. +GRID 85 .7 .7 0. +GRID 86 .7 .8 0. +GRID 87 .7 .9 0. +GRID 88 .7 1. 0. +GRID 89 .8 0. 0. +GRID 90 .8 .1 0. +GRID 91 .8 .2 0. +GRID 92 .8 .3 0. +GRID 93 .8 .4 0. +GRID 94 .8 .5 0. +GRID 95 .8 .6 0. +GRID 96 .8 .7 0. +GRID 97 .8 .8 0. +GRID 98 .8 .9 0. +GRID 99 .8 1. 0. +GRID 100 .9 0. 0. +GRID 101 .9 .1 0. +GRID 102 .9 .2 0. +GRID 103 .9 .3 0. +GRID 104 .9 .4 0. +GRID 105 .9 .5 0. +GRID 106 .9 .6 0. +GRID 107 .9 .7 0. +GRID 108 .9 .8 0. +GRID 109 .9 .9 0. +GRID 110 .9 1. 0. +GRID 111 1. 0. 0. +GRID 112 1. .1 0. +GRID 113 1. .2 0. +GRID 114 1. .3 0. +GRID 115 1. .4 0. +GRID 116 1. .5 0. +GRID 117 1. .6 0. +GRID 118 1. .7 0. +GRID 119 1. .8 0. +GRID 120 1. .9 0. +GRID 121 1. 1. 0. +$ELEMENTS +CQUAD4 1 1 1 12 13 2 +CQUAD4 2 2 2 13 14 3 +CQUAD4 3 3 3 14 15 4 +CQUAD4 4 4 4 15 16 5 +CQUAD4 5 5 5 16 17 6 +CQUAD4 6 6 6 17 18 7 +CQUAD4 7 7 7 18 19 8 +CQUAD4 8 8 8 19 20 9 +CQUAD4 9 9 9 20 21 10 +CQUAD4 10 10 10 21 22 11 +CQUAD4 11 11 12 23 24 13 +CQUAD4 12 12 13 24 25 14 +CQUAD4 13 13 14 25 26 15 +CQUAD4 14 14 15 26 27 16 +CQUAD4 15 15 16 27 28 17 +CQUAD4 16 16 17 28 29 18 +CQUAD4 17 17 18 29 30 19 +CQUAD4 18 18 19 30 31 20 +CQUAD4 19 19 20 31 32 21 +CQUAD4 20 20 21 32 33 22 +CQUAD4 21 21 23 34 35 24 +CQUAD4 22 22 24 35 36 25 +CQUAD4 23 23 25 36 37 26 +CQUAD4 24 24 26 37 38 27 +CQUAD4 25 25 27 38 39 28 +CQUAD4 26 26 28 39 40 29 +CQUAD4 27 27 29 40 41 30 +CQUAD4 28 28 30 41 42 31 +CQUAD4 29 29 31 42 43 32 +CQUAD4 30 30 32 43 44 33 +CQUAD4 31 31 34 45 46 35 +CQUAD4 32 32 35 46 47 36 +CQUAD4 33 33 36 47 48 37 +CQUAD4 34 34 37 48 49 38 +CQUAD4 35 35 38 49 50 39 +CQUAD4 36 36 39 50 51 40 +CQUAD4 37 37 40 51 52 41 +CQUAD4 38 38 41 52 53 42 +CQUAD4 39 39 42 53 54 43 +CQUAD4 40 40 43 54 55 44 +CQUAD4 41 41 45 56 57 46 +CQUAD4 42 42 46 57 58 47 +CQUAD4 43 43 47 58 59 48 +CQUAD4 44 44 48 59 60 49 +CQUAD4 45 45 49 60 61 50 +CQUAD4 46 46 50 61 62 51 +CQUAD4 47 47 51 62 63 52 +CQUAD4 48 48 52 63 64 53 +CQUAD4 49 49 53 64 65 54 +CQUAD4 50 50 54 65 66 55 +CQUAD4 51 51 56 67 68 57 +CQUAD4 52 52 57 68 69 58 +CQUAD4 53 53 58 69 70 59 +CQUAD4 54 54 59 70 71 60 +CQUAD4 55 55 60 71 72 61 +CQUAD4 56 56 61 72 73 62 +CQUAD4 57 57 62 73 74 63 +CQUAD4 58 58 63 74 75 64 +CQUAD4 59 59 64 75 76 65 +CQUAD4 60 60 65 76 77 66 +CQUAD4 61 61 67 78 79 68 +CQUAD4 62 62 68 79 80 69 +CQUAD4 63 63 69 80 81 70 +CQUAD4 64 64 70 81 82 71 +CQUAD4 65 65 71 82 83 72 +CQUAD4 66 66 72 83 84 73 +CQUAD4 67 67 73 84 85 74 +CQUAD4 68 68 74 85 86 75 +CQUAD4 69 69 75 86 87 76 +CQUAD4 70 70 76 87 88 77 +CQUAD4 71 71 78 89 90 79 +CQUAD4 72 72 79 90 91 80 +CQUAD4 73 73 80 91 92 81 +CQUAD4 74 74 81 92 93 82 +CQUAD4 75 75 82 93 94 83 +CQUAD4 76 76 83 94 95 84 +CQUAD4 77 77 84 95 96 85 +CQUAD4 78 78 85 96 97 86 +CQUAD4 79 79 86 97 98 87 +CQUAD4 80 80 87 98 99 88 +CQUAD4 81 81 89 100 101 90 +CQUAD4 82 82 90 101 102 91 +CQUAD4 83 83 91 102 103 92 +CQUAD4 84 84 92 103 104 93 +CQUAD4 85 85 93 104 105 94 +CQUAD4 86 86 94 105 106 95 +CQUAD4 87 87 95 106 107 96 +CQUAD4 88 88 96 107 108 97 +CQUAD4 89 89 97 108 109 98 +CQUAD4 90 90 98 109 110 99 +CQUAD4 91 91 100 111 112 101 +CQUAD4 92 92 101 112 113 102 +CQUAD4 93 93 102 113 114 103 +CQUAD4 94 94 103 114 115 104 +CQUAD4 95 95 104 115 116 105 +CQUAD4 96 96 105 116 117 106 +CQUAD4 97 97 106 117 118 107 +CQUAD4 98 98 107 118 119 108 +CQUAD4 99 99 108 119 120 109 +CQUAD4 100 100 109 120 121 110 +$SPCs +SPC1 1 123456 1 11 +SPC1 1 123456 12 22 +SPC1 1 123456 23 33 +SPC1 1 123456 34 44 +SPC1 1 123456 45 55 +SPC1 1 123456 56 66 +SPC1 1 123456 67 77 +SPC1 1 123456 78 88 +SPC1 1 123456 89 99 +SPC1 1 123456 100 110 +SPC1 1 123456 111 121 +SPC1 1 123456 1 111 +SPC1 1 123456 2 112 +SPC1 1 123456 3 113 +SPC1 1 123456 4 114 +SPC1 1 123456 5 115 +SPC1 1 123456 6 116 +SPC1 1 123456 7 117 +SPC1 1 123456 8 118 +SPC1 1 123456 9 119 +SPC1 1 123456 10 120 +SPC1 1 123456 11 121