diff --git a/CODEOWNERS b/CODEOWNERS
new file mode 100644
index 0000000000000000000000000000000000000000..90d402033f6ec302a2ca21305c37cf8366a628be
--- /dev/null
+++ b/CODEOWNERS
@@ -0,0 +1,16 @@
+# This is a comment.
+# Each line is a file pattern followed by one or more owners.
+
+# These owners will be the default owners for everything in
+# the repo. Unless a later match takes precedence, they
+# will be requested for review when someone opens a pull request.
+*       @R.Boman
+*       @acrovato
+
+# These owners own a directory and nested subdirectories
+
+# You can also use email addresses if you prefer. They'll be
+# used to look up users just like we do for commit author
+# emails.
+#*.go docs@example.com
+
diff --git a/CREDITS.txt b/CREDITS.txt
new file mode 100644
index 0000000000000000000000000000000000000000..eceb05fadb44bc4a84d3d03b6ff4eefb8c21fce4
--- /dev/null
+++ b/CREDITS.txt
@@ -0,0 +1,3 @@
+dartflo_mphys is Copyright (C) of University of Liège.
+
+The main author of dartflo_mphys is Adrien Crovato.
diff --git a/LICENSE.txt b/LICENSE.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8dada3edaf50dbc082c9a125058f25def75e625a
--- /dev/null
+++ b/LICENSE.txt
@@ -0,0 +1,201 @@
+                                 Apache License
+                           Version 2.0, January 2004
+                        http://www.apache.org/licenses/
+
+   TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
+
+   1. Definitions.
+
+      "License" shall mean the terms and conditions for use, reproduction,
+      and distribution as defined by Sections 1 through 9 of this document.
+
+      "Licensor" shall mean the copyright owner or entity authorized by
+      the copyright owner that is granting the License.
+
+      "Legal Entity" shall mean the union of the acting entity and all
+      other entities that control, are controlled by, or are under common
+      control with that entity. For the purposes of this definition,
+      "control" means (i) the power, direct or indirect, to cause the
+      direction or management of such entity, whether by contract or
+      otherwise, or (ii) ownership of fifty percent (50%) or more of the
+      outstanding shares, or (iii) beneficial ownership of such entity.
+
+      "You" (or "Your") shall mean an individual or Legal Entity
+      exercising permissions granted by this License.
+
+      "Source" form shall mean the preferred form for making modifications,
+      including but not limited to software source code, documentation
+      source, and configuration files.
+
+      "Object" form shall mean any form resulting from mechanical
+      transformation or translation of a Source form, including but
+      not limited to compiled object code, generated documentation,
+      and conversions to other media types.
+
+      "Work" shall mean the work of authorship, whether in Source or
+      Object form, made available under the License, as indicated by a
+      copyright notice that is included in or attached to the work
+      (an example is provided in the Appendix below).
+
+      "Derivative Works" shall mean any work, whether in Source or Object
+      form, that is based on (or derived from) the Work and for which the
+      editorial revisions, annotations, elaborations, or other modifications
+      represent, as a whole, an original work of authorship. For the purposes
+      of this License, Derivative Works shall not include works that remain
+      separable from, or merely link (or bind by name) to the interfaces of,
+      the Work and Derivative Works thereof.
+
+      "Contribution" shall mean any work of authorship, including
+      the original version of the Work and any modifications or additions
+      to that Work or Derivative Works thereof, that is intentionally
+      submitted to Licensor for inclusion in the Work by the copyright owner
+      or by an individual or Legal Entity authorized to submit on behalf of
+      the copyright owner. For the purposes of this definition, "submitted"
+      means any form of electronic, verbal, or written communication sent
+      to the Licensor or its representatives, including but not limited to
+      communication on electronic mailing lists, source code control systems,
+      and issue tracking systems that are managed by, or on behalf of, the
+      Licensor for the purpose of discussing and improving the Work, but
+      excluding communication that is conspicuously marked or otherwise
+      designated in writing by the copyright owner as "Not a Contribution."
+
+      "Contributor" shall mean Licensor and any individual or Legal Entity
+      on behalf of whom a Contribution has been received by Licensor and
+      subsequently incorporated within the Work.
+
+   2. Grant of Copyright License. Subject to the terms and conditions of
+      this License, each Contributor hereby grants to You a perpetual,
+      worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+      copyright license to reproduce, prepare Derivative Works of,
+      publicly display, publicly perform, sublicense, and distribute the
+      Work and such Derivative Works in Source or Object form.
+
+   3. Grant of Patent License. Subject to the terms and conditions of
+      this License, each Contributor hereby grants to You a perpetual,
+      worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+      (except as stated in this section) patent license to make, have made,
+      use, offer to sell, sell, import, and otherwise transfer the Work,
+      where such license applies only to those patent claims licensable
+      by such Contributor that are necessarily infringed by their
+      Contribution(s) alone or by combination of their Contribution(s)
+      with the Work to which such Contribution(s) was submitted. If You
+      institute patent litigation against any entity (including a
+      cross-claim or counterclaim in a lawsuit) alleging that the Work
+      or a Contribution incorporated within the Work constitutes direct
+      or contributory patent infringement, then any patent licenses
+      granted to You under this License for that Work shall terminate
+      as of the date such litigation is filed.
+
+   4. Redistribution. You may reproduce and distribute copies of the
+      Work or Derivative Works thereof in any medium, with or without
+      modifications, and in Source or Object form, provided that You
+      meet the following conditions:
+
+      (a) You must give any other recipients of the Work or
+          Derivative Works a copy of this License; and
+
+      (b) You must cause any modified files to carry prominent notices
+          stating that You changed the files; and
+
+      (c) You must retain, in the Source form of any Derivative Works
+          that You distribute, all copyright, patent, trademark, and
+          attribution notices from the Source form of the Work,
+          excluding those notices that do not pertain to any part of
+          the Derivative Works; and
+
+      (d) If the Work includes a "NOTICE" text file as part of its
+          distribution, then any Derivative Works that You distribute must
+          include a readable copy of the attribution notices contained
+          within such NOTICE file, excluding those notices that do not
+          pertain to any part of the Derivative Works, in at least one
+          of the following places: within a NOTICE text file distributed
+          as part of the Derivative Works; within the Source form or
+          documentation, if provided along with the Derivative Works; or,
+          within a display generated by the Derivative Works, if and
+          wherever such third-party notices normally appear. The contents
+          of the NOTICE file are for informational purposes only and
+          do not modify the License. You may add Your own attribution
+          notices within Derivative Works that You distribute, alongside
+          or as an addendum to the NOTICE text from the Work, provided
+          that such additional attribution notices cannot be construed
+          as modifying the License.
+
+      You may add Your own copyright statement to Your modifications and
+      may provide additional or different license terms and conditions
+      for use, reproduction, or distribution of Your modifications, or
+      for any such Derivative Works as a whole, provided Your use,
+      reproduction, and distribution of the Work otherwise complies with
+      the conditions stated in this License.
+
+   5. Submission of Contributions. Unless You explicitly state otherwise,
+      any Contribution intentionally submitted for inclusion in the Work
+      by You to the Licensor shall be under the terms and conditions of
+      this License, without any additional terms or conditions.
+      Notwithstanding the above, nothing herein shall supersede or modify
+      the terms of any separate license agreement you may have executed
+      with Licensor regarding such Contributions.
+
+   6. Trademarks. This License does not grant permission to use the trade
+      names, trademarks, service marks, or product names of the Licensor,
+      except as required for reasonable and customary use in describing the
+      origin of the Work and reproducing the content of the NOTICE file.
+
+   7. Disclaimer of Warranty. Unless required by applicable law or
+      agreed to in writing, Licensor provides the Work (and each
+      Contributor provides its Contributions) on an "AS IS" BASIS,
+      WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
+      implied, including, without limitation, any warranties or conditions
+      of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A
+      PARTICULAR PURPOSE. You are solely responsible for determining the
+      appropriateness of using or redistributing the Work and assume any
+      risks associated with Your exercise of permissions under this License.
+
+   8. Limitation of Liability. In no event and under no legal theory,
+      whether in tort (including negligence), contract, or otherwise,
+      unless required by applicable law (such as deliberate and grossly
+      negligent acts) or agreed to in writing, shall any Contributor be
+      liable to You for damages, including any direct, indirect, special,
+      incidental, or consequential damages of any character arising as a
+      result of this License or out of the use or inability to use the
+      Work (including but not limited to damages for loss of goodwill,
+      work stoppage, computer failure or malfunction, or any and all
+      other commercial damages or losses), even if such Contributor
+      has been advised of the possibility of such damages.
+
+   9. Accepting Warranty or Additional Liability. While redistributing
+      the Work or Derivative Works thereof, You may choose to offer,
+      and charge a fee for, acceptance of support, warranty, indemnity,
+      or other liability obligations and/or rights consistent with this
+      License. However, in accepting such obligations, You may act only
+      on Your own behalf and on Your sole responsibility, not on behalf
+      of any other Contributor, and only if You agree to indemnify,
+      defend, and hold each Contributor harmless for any liability
+      incurred by, or claims asserted against, such Contributor by reason
+      of your accepting any such warranty or additional liability.
+
+   END OF TERMS AND CONDITIONS
+
+   APPENDIX: How to apply the Apache License to your work.
+
+      To apply the Apache License to your work, attach the following
+      boilerplate notice, with the fields enclosed by brackets "{}"
+      replaced with your own identifying information. (Don't include
+      the brackets!)  The text should be enclosed in the appropriate
+      comment syntax for the file format. We also recommend that a
+      file or class name and description of purpose be included on the
+      same "printed page" as the copyright notice for easier
+      identification within third-party archives.
+
+   Copyright {yyyy} {name of copyright owner}
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.
diff --git a/README.md b/README.md
index 7d8e07e105d5cc6c15f775b3ba09f1ad013d5f0a..8666a4a90e98cab8f8faadf8f363e007b15940ac 100644
--- a/README.md
+++ b/README.md
@@ -1,5 +1,5 @@
 # DARTFlo - MPhys
-This repository contains some examples of analysis and optimization cases with [DART](https://gitlab.uliege.be/am-dept/dartflo) and [OpenMDAO](https://openmdao.org/)/[MPhys](https://github.com/OpenMDAO/mphys).
+This repository contains some examples of analysis and optimization cases with [DART](https://gitlab.uliege.be/am-dept/dartflo) and [OpenMDAO](https://openmdao.org/)/[MPhys](https://github.com/openmdao/mphys).
 
 ## List of cases
 - crm: aerostructural analysis of the NASA CRM.
@@ -11,13 +11,13 @@ This repository contains some examples of analysis and optimization cases with [
 ## Version of solvers
 - [OpenMDAO](https://openmdao.org/): v3.30
 - [pyOptSparse](https://github.com/mdolab/pyoptsparse): v2.10.1
-- [MPhys](https://github.com/OpenMDAO/mphys): v1.3.0
+- [MPhys](https://github.com/openmdao/mphys): v2.0.0
 - [OMFlut](https://gitlab.uliege.be/am-dept/omflut): v1.1.0
 - [pySpline](https://github.com/mdolab/pyspline): v1.5.1
 - [pyGeo](https://github.com/mdolab/pygeo): v1.13.0
 - [Gmsh](https://gmsh.info/): v4.10.5
-- [DART](https://gitlab.uliege.be/am-dept/dartflo): v1.2.0
-- [MELD](https://github.com/smdogroup/funtofem): git revision f2b39ef (06/28/2022)
-- [TACS](https://github.com/smdogroup/tacs): v3.7.1
+- [DART](https://gitlab.uliege.be/am-dept/dartflo): v1.3.0
+- [MELD](https://github.com/smdogroup/funtofem): v0.3.9
+- [TACS](https://github.com/smdogroup/tacs): v3.9.1
 - [SDPM](https://gitlab.uliege.be/am-dept/sdpm): v1.2.0
 - [PyPk](https://gitlab.uliege.be/am-dept/pypk): v1.1.0
diff --git a/crm/geometry.py b/crm/geometry.py
index a0bd575896d727ffe5d1c1e444fbd6de27f99924..caf055557a6fe1d5b849dcc3cbedfb8eec9fbac2 100644
--- a/crm/geometry.py
+++ b/crm/geometry.py
@@ -4,7 +4,13 @@
 import numpy as np
 from mpi4py import MPI
 import openmdao.api as om
-from mphys.builder import Builder
+from mphys import Builder, MPhysVariables
+
+# Define short aliases for convenience
+XAERO_IN = MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_INPUT
+XAERO_OUT = MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_OUTPUT
+XSTRUCT_IN = MPhysVariables.Structures.Geometry.COORDINATES_INPUT
+XSTRUCT_OUT = MPhysVariables.Structures.Geometry.COORDINATES_OUTPUT
 
 # Geometry parametrization and constraints
 class PyGeoGeom(om.ExplicitComponent):
@@ -20,9 +26,12 @@ class PyGeoGeom(om.ExplicitComponent):
         self.omPtSetList = []
         # Setup geometry, constraints and their I/O
         # add IO for surface meshes, the DVGeo object must be configured later by calling addDisciplinePoints
-        for disc in self.config['disciplines']:
-            self.add_input(f'x_{disc}_in', distributed=True, shape_by_conn=True)
-            self.add_output(f'x_{disc}0', distributed=True, copy_shape=f'x_{disc}_in', tags=['mphys_coordinates'])
+        if 'aero' in self.config['disciplines']:
+            self.add_input(XAERO_IN, distributed=True, shape_by_conn=True, tags=['mphys_coordinates'])
+            self.add_output(XAERO_OUT, distributed=True, copy_shape=XAERO_IN, tags=['mphys_coordinates'])
+        if 'struct' in self.config['disciplines']:
+            self.add_input(XSTRUCT_IN, distributed=True, shape_by_conn=True, tags=['mphys_coordinates'])
+            self.add_output(XSTRUCT_OUT, distributed=True, copy_shape=XSTRUCT_IN, tags=['mphys_coordinates'])
         # add a reference axis
         if 'refAxis' in self.config:
             for name, prm in self.config['refAxis'].items():
@@ -88,8 +97,12 @@ class PyGeoGeom(om.ExplicitComponent):
     def addDisciplinePoints(self, discipline, points, **kwargs):
         """Add the points of the discipline to be embedded in the parametrization
         """
-        self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), f'x_{discipline}0', **kwargs)
-        self.omPtSetList.append(f'x_{discipline}0')
+        if discipline == 'aero':
+            self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), XAERO_OUT, **kwargs)
+            self.omPtSetList.append(XAERO_OUT)
+        elif discipline == 'struct':
+            self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), XSTRUCT_OUT, **kwargs)
+            self.omPtSetList.append(XSTRUCT_OUT)
 
     def setConstrainedSurface(self, surface):
         """Set (a triangulation of) the surface to be constrained and add the constraints to DVCon
diff --git a/crm/run.py b/crm/run.py
index 9b36b7037347481ba1067b389029f802226d6060..6b41286b2d979d82aa95647c77c7f401ccb818b0 100644
--- a/crm/run.py
+++ b/crm/run.py
@@ -147,7 +147,7 @@ def main():
             matrix[i,7] = case.get_constraints(scaled=False)['flight.trim_cruise.n'][0]
             matrix[i,8] = case.get_constraints(scaled=False)['flight.trim_maneuver.n'][0]
             matrix[i,9] = case.get_constraints(scaled=False)['static.maneuver.ks_fail'][0]
-            matrix[i,10] = case.get_constraints(scaled=False)['flight.volume_match.vc'][0]
+            #matrix[i,10] = case.get_constraints(scaled=False)['flight.volume_match.vc'][0]
             matrix[i,11] = case.get_constraints(scaled=False)['flight.fuel_match.fm'][0]
             matrix[i,12] = case.get_objectives(scaled=False)['flight.fuel_burn.fb'][0]
         np.savetxt('history.dat', matrix, header='it, aoa0, aoa1, cl0, cl1, cd0, ww, n0, n1, ks1, vol, fm, wf')
diff --git a/crm/setup/crm.geo b/crm/setup/crm.geo
index e0085f9366f93ed47f11bb1336612ae4f6838e7d..270847abe2b5f93d2182eefced1222dfc1dab55b 100644
--- a/crm/setup/crm.geo
+++ b/crm/setup/crm.geo
@@ -1,123 +1,167 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-/* NASA CRM - NASA Common Research Model
- * CAD created by Guillem Battle I Capa
- * Valid mesh obtained with Gmsh 4.10.5
- */
-
-/* Model information
- * Length of fuselage: 62.79 m
- * Chord length at wing root: 12.09 m
- * Chord length at Yehudi break: 7.30 m
- * Chord length at wing tip: 2.76 m
- * Chord length at tail root: 5.45 m
- * Chord length at tail tip: 2.33 m
- * Reference wing area: 1/2 * 383.69 m^2
- * Reference chord length: 7.01 m
- */
-
-/// Constants
-scl = 2; // scaling factor for lifting surfaces mesh size (testing purpose)
-DefineConstant[ msFar = {36.8, Name "Mesh size on farfield" },
-                msFus = {0.628, Name "Mesh size on fuselage" },
-                msWingR = {0.121 * scl, Name "Mesh size on wing root" },
-				msWingY = {0.073 * scl, Name "Mesh size on wing Yehudi break" },
-				msWingT = {0.028 * scl, Name "Mesh size on wing tip" },
-				msTailR = {0.055 * scl, Name "Mesh size on tail root" },
-                msTailT = {0.023 * scl, Name "Mesh size on tail tip" } ];
-
-/// Geometry definition
-// Import
-Geometry.OCCAutoFix = 1;
-Geometry.OCCFixDegenerated = 1;
-Geometry.OCCFixSmallEdges = 1;
-Geometry.OCCFixSmallFaces = 1;
-Geometry.OCCSewFaces = 1;
-Geometry.Tolerance = 1e-5; // 0.01 mm
-Geometry.OCCTargetUnit = "M";
-Merge "crm.stp";
-Coherence;
-
-// Re-orient
-ReverseMesh Surface{1:48}; //To reverse the direction of the mesh normal vectors
-
-// Volume
-Surface Loop(1) = {4:48};
-Volume(1) = {1};
-
-// Embedded
-Line{3} In Surface{16};
-Line{7} In Surface{16};
-Line{15} In Surface{16};
-Surface{1} In Volume{1};
-Surface{2} In Volume{1};
-Surface{3} In Volume{1};
-
-// Physical groups
-Physical Line("wakeExW") = {10,11,12,13};
-Physical Line("teW") = {9,14};
-Physical Line("wakeTipW") = {10};
-Physical Line("wakeExT") = {1,2,5};
-Physical Line("teT") = {6};
-Physical Line("wakeTipT") = {1,2};
-Physical Surface("symmetry") = {19,20,21,22};
-Physical Surface("downstream") = {16};
-Physical Surface("upstream") = {18};
-Physical Surface("farfield") = {14,15,17};
-Physical Surface("fuselage") = {13,23:45};
-Physical Surface("wing") = {4,5,6,11,12};
-Physical Surface("wing_") = {7,8,9,10};
-Physical Surface("tail") = {47,48};
-Physical Surface("tail_") = {46};
-Physical Surface("wakeW") = {2,3};
-Physical Surface("wakeT") = {1};
-Physical Volume("field") = {1};
-
-/// Mesh definition
-// Farfield
-Characteristic Length {34,35,36,37,38,39,40,43,44,45} = msFar; // farfield boundaries
-Characteristic Length {3,4,7,8,15,41,42} = msFar; // wake
-
-// Fuselage
-Characteristic Length {11,12,13,51,61,50,73,74,62,63,64,65,47,29,70,46} = msFus; // body
-Characteristic Length {51,52,53,54,55,56,59,60,89,85,86,90,91,92} = 0.5 * msFus; // nose and windshield
-Characteristic Length {57,58,87,88} = 0.2 * msFus; // surface between windshield
-Characteristic Length {49,72,28,33,32,24,31,30,75,71,48,76,69,47} = 0.75 * msFus; // wingbox 
-Characteristic Length {5,66,67,68} = 0.25 * msFus; // APU
-Characteristic Length {46,11} = msFus; // Fuselage-wake intersection
-
-// Wing
-Characteristic Length {24, 25,26,27,14} = msWingR; // wing root
-Characteristic Length {22, 20,21,23,9} = msWingY; // wing kink
-Characteristic Length {18,10,16,19} = 3 * msWingT; // wing tip TE
-Characteristic Length {17} = msWingT; // wing tip LE
-
-// Tail
-Characteristic Length {79:84,93:99, 100:102,77,78,103,6} = msTailR; // Tail root
-Characteristic Length {1,2,105} = 2 * msTailT; // tail tip TE
-Characteristic Length {104} = msTailT; // tail tip LE
-
-// Mesh algos
-Mesh.Algorithm = 6;
-Mesh.Algorithm3D = 10;
-Mesh.Optimize = 1;
-Mesh.Smoothing = 10;
-Mesh.SmoothNormals = 1;
-
-Mesh 1;
-Geometry.Tolerance = 1e-6;
-Coherence Mesh;
+// NASA CRM - NASA Common Research Model
+// CAD created by Guillem Battle I Capa and Adrien Crovato
+// Valid mesh obtained with gmsh 4.10.5
+
+// Model information
+// Length of fuselage: 62.79 m
+// Chord length at wing root: 12.09 m
+// Chord length at Yehudi break: 7.30 m
+// Chord length at wing tip: 2.76 m
+// Chord length at tail root: 5.45 m
+// Chord length at tail tip: 2.33 m
+// Reference wing area: 1/2 * 383.69 m^2
+// Reference chord length: 7.01 m
+
+/// Constants
+DefineConstant[ surfMeshSize = {0.01, Name "Factor multiplying the chord length to obtain mesh size on the body surface" },
+				growthRatio = {1.2, Name "Growth ratio for elements in the field" }
+				multTipTe = {5, Name "Mesh size multiplicator for wing tip trailing edge" } ];
+
+// OCC options
+SetFactory("OpenCASCADE");
+Geometry.OCCAutoFix = 1;
+Geometry.OCCFixDegenerated = 1;
+Geometry.OCCFixSmallEdges = 1;
+Geometry.OCCFixSmallFaces = 1;
+Geometry.OCCSewFaces = 1;
+Geometry.Tolerance = 1e-5;
+Geometry.OCCTargetUnit = "M";
+
+// CRM geometry
+Merge "crm.stp";
+Geometry.Tolerance = 1e-6;
+
+// Add wing wake
+Point(1001) = {1260,     0, 3.10};
+Point(1002) = {1260,  2.14, 3.38};
+Point(1003) = {1260,  3.07, 3.56};
+Point(1004) = {1260, 10.87, 4.49};
+Point(1005) = {1260, 29.38, 6.87};
+Line(1001) = {53, 1001};
+Line(1002) = {37, 1005};
+Line(2001) = {1001, 1002};
+Line(2002) = {1002, 1003};
+Line(2003) = {1003, 1004};
+Line(2004) = {1004, 1005};
+Curve Loop(1001) = {1001, 2001, 2002, 2003, 2004, 1002, 53, 62, 69, 80, 79};
+Surface(37) = {1001};
+
+// Add tail wake
+Point(1011) = {1260,     0, 6.49};
+Point(1012) = {1260,  0.85, 6.57};
+Point(1013) = {1260, 10.65, 7.67};
+Line(1011) = {64, 1011};
+Line(1012) = {78, 1013};
+Line(2011) = {1011, 1012};
+Line(2012) = {1012, 1013};
+Curve Loop(1011) = {1011, 2011, 2012, 1012, 111, 92};
+Surface(38) = {1011};
+
+// Domain
+Point(100001) = {-1260, 0., -1260};
+Point(100002) = {1260, 0., -1260};
+Point(100003) = {1260, 0., 1260};
+Point(100004) = {-1260, 0., 1260};
+Point(100011) = {-1260, 1260, -1260};
+Point(100012) = {1260, 1260, -1260};
+Point(100013) = {1260, 1260, 1260};
+Point(100014) = {-1260, 1260, 1260};
+Line(100001) = {100001, 100002};
+Line(100002) = {100002, 1001};
+Line(100003) = {1001, 1011};
+Line(100004) = {1011, 100003};
+Line(100005) = {100003, 100004};
+Line(100006) = {100004, 100001};
+Line(100011) = {100011, 100012};
+Line(100012) = {100012, 100013};
+Line(100013) = {100013, 100014};
+Line(100014) = {100014, 100011};
+Line(100021) = {100001, 100011};
+Line(100022) = {100002, 100012};
+Line(100023) = {100003, 100013};
+Line(100024) = {100004, 100014};
+Curve Loop(100001) = {100021, 100014, 100024, 100006};
+Curve Loop(100002) = {100022, 100012, 100023, 100004, 100003, 100002};
+Curve Loop(100003) = {100001, 100022, 100011, 100021};
+Curve Loop(100004) = {100005, 100024, 100013, 100023};
+Curve Loop(100005) = {100011, 100012, 100013, 100014};
+Curve Loop(100006) = {100001, 100002, 100003, 100004, 100005, 100006};
+Curve Loop(100007) = {16, 9, 7, 3, 22, 25, 93, 71, 46, 41, 90, 91, 108, 109, 77, 76, 30, 35, 94, 26, 23, 19, 15};
+Plane Surface(100001) = {100001};
+Plane Surface(100002) = {100002};
+Plane Surface(100003) = {100003};
+Plane Surface(100004) = {100004};
+Plane Surface(100005) = {100005};
+Plane Surface(100006) = {100006, 100007};
+Surface Loop(1) = {1:36, 100001:100006};
+Volume(1) = {1};
+
+// Embedded
+Curve{2001} In Surface{100002};
+Curve{2002} In Surface{100002};
+Curve{2003} In Surface{100002};
+Curve{2004} In Surface{100002};
+Curve{2011} In Surface{100002};
+Curve{2012} In Surface{100002};
+Curve{1001} In Surface{100006};
+Curve{1011} In Surface{100006};
+Surface{37} In Volume{1};
+Surface{38} In Volume{1};
+
+// Physical groups
+Physical Line("wingTe") = {53, 62};
+Physical Line("wingWakeTip") = {1002};
+Physical Line("wingWakeEx") = {69, 80, 79};
+Physical Line("tailTe") = {111};
+Physical Line("tailWakeTip") = {1012};
+Physical Line("tailWakeEx") = {92};
+Physical Surface("symmetry") = {100006};
+Physical Surface("downstream") = {100002};
+Physical Surface("upstream") = {100001};
+Physical Surface("farfield") = {100003, 100004, 100005};
+Physical Surface("fuselage") = {1:14, 24:33};
+Physical Surface("wing") = {22, 23, 16, 17, 15};
+Physical Surface("wing_") = {20, 21, 18, 19};
+Physical Surface("tail") = {35, 36};
+Physical Surface("tail_") = {34};
+Physical Surface("wingWake") = {37};
+Physical Surface("tailWake") = {38};
+Physical Volume("field") = {1};
+
+// Mesh sizes
+msFus = 62.8 * surfMeshSize;
+msWingR = 12.1 * surfMeshSize;
+msWingY = 7.3 * surfMeshSize;
+msWingT = 2.8 * surfMeshSize;
+msTailR = 5.5 * surfMeshSize;
+msTailT = 2.3 * surfMeshSize;
+msFar = msWingR * growthRatio ^ (Log(1 - (1 - growthRatio) * 1260 / msWingR) / Log(growthRatio) - 1);
+
+Characteristic Length {100001:100006, 100011:100014} = msFar; // farfield
+Characteristic Length {1001:1005, 1011:1013} = 30; // wake
+
+Characteristic Length {16:19, 49, 28, 25, 26, 34, 22, 30, 23, 32, 33, 54, 51, 50, 53} = msFus; // body
+Characteristic Length {1:3, 8:15} = 0.5 * msFus; // nose and windshield
+Characteristic Length {4:7} = 0.2 * msFus; // surface between windshield
+Characteristic Length {20, 24, 29, 27, 31, 47, 48, 52, 21} = 0.75 * msFus; // wingbox
+Characteristic Length {63, 64, 76} = 0.25 * msFus; // APU
+
+Characteristic Length {43:46} = msWingR; // wing root
+Characteristic Length {39:42} = msWingY; // wing kink
+Characteristic Length {36} = msWingT; // wing tip LE
+Characteristic Length {35, 37, 38} = multTipTe * msWingT; // wing tip TE
+
+Characteristic Length {55:62, 65:75} = msTailR; // Tail root
+Characteristic Length {77} = msTailT; // tail tip LE
+Characteristic Length {78} = 2 * msTailT; // tail tip TE
+
+// Mesh control
+Transfinite Curve{48, 49} = 2;
+ReverseMesh Surface{1:36};
+ReverseMesh Surface{100002, 100006};
+Mesh.Algorithm = 6;
+Mesh.Algorithm3D = 10;
+Mesh.OptimizeNetgen = 0;
+Mesh.Optimize = 1;
+Mesh.Smoothing = 10;
+Mesh.SmoothNormals = 1;
diff --git a/crm/setup/crm.stp b/crm/setup/crm.stp
index b50b1bb88fde5d9db7ac23ff5b4a6605285eed8b..f498df840d7c7feb58e757c1a65244d40a29850d 100644
--- a/crm/setup/crm.stp
+++ b/crm/setup/crm.stp
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:4156328b8ca196e81de869b1144d7b21161e07e5a8a4bd4d89f1a3a447f4bca2
-size 16233634
+oid sha256:aceebaa286b096fff2e2e93776ae0509a34ad37dfb53a61398fff577dfd4bcba
+size 20941895
diff --git a/crm/setup/dart_setup.py b/crm/setup/dart_setup.py
index 1bd4894b7086767879e7de5a9f0e3631509a7e56..e590dce344542b51234e08461b5476ab570ff92a 100644
--- a/crm/setup/dart_setup.py
+++ b/crm/setup/dart_setup.py
@@ -8,9 +8,7 @@ def get(nthrd, verb, iscn):
         'Verb' : verb, # verbosity
         # Model (geometry or mesh)
         'File' : get_path('crm.geo', 'setup'), # Input file containing the model
-        'Pars' : {'msWingR': 0.121, 'msWingY': 0.073, 'msWingT': 0.028,
-                  'msTailR': 0.055, 'msTailT': 0.023,
-                  'msFus': 0.628, 'msFar': 36.8}, # parameters for input file model
+        'Pars' : {'surfMeshSize': 0.01, 'growthRatio': 1.2, 'multTipTe': 5.5}, # parameters for input file model
         'Dim' : 3, # problem dimension
         'Format' : 'vtk', # save format (vtk or gmsh)
         # Markers
@@ -19,10 +17,10 @@ def get(nthrd, verb, iscn):
         'Fuselage' : 'fuselage', # Name of physical group containing the fuselage boundaries
         'Symmetry' : 'symmetry',
         'Wings' : ['wing', 'tail'], # LIST of names of physical groups containing the lifting surface body
-        'Wakes' : ['wakeW', 'wakeT'], # LIST of names of physical group containing the wake
-        'WakeTips' : ['wakeTipW', 'wakeTipT'], # LIST of names of physical group containing the edge of the wake
-        'WakeExs' : ['wakeExW', 'wakeExT'], # LIST of names of physical group containing the edge of the wake and the intersection of lifting surface with fuselage (to be excluded from Wake B.C.)
-        'Tes' : ['teW', 'teT'], # LIST of names of physical group containing the trailing edge
+        'Wakes' : ['wingWake', 'tailWake'], # LIST of names of physical group containing the wake
+        'WakeTips' : ['wingWakeTip', 'tailWakeTip'], # LIST of names of physical group containing the edge of the wake
+        'WakeExs' : ['wingWakeEx', 'tailWakeEx'], # LIST of names of physical group containing the edge of the wake and the intersection of lifting surface with fuselage (to be excluded from Wake B.C.)
+        'Tes' : ['wingTe', 'tailTe'], # LIST of names of physical group containing the trailing edge
         # Freestream
         #'M_inf' : 0.85, # freestream Mach number
         #'AoA' : 1.0, # freestream angle of attack
diff --git a/crm/setup/meld_setup.py b/crm/setup/meld_setup.py
new file mode 100644
index 0000000000000000000000000000000000000000..ea3ed9eb08bc384511e18abf3082ad5f4cef5673
--- /dev/null
+++ b/crm/setup/meld_setup.py
@@ -0,0 +1,7 @@
+def get():
+    return {
+        'body_tags': [{
+            'aero': None,
+            'struct': ['Skin', 'Spar']
+            }]
+    }
diff --git a/crm/setup/tacs_setup.py b/crm/setup/tacs_setup.py
index fec793bb51dd7b687bd29047288bf0717518c392..01f7cb3bf6edc68b9b19c1dc6afe01178372ec7a 100644
--- a/crm/setup/tacs_setup.py
+++ b/crm/setup/tacs_setup.py
@@ -1,5 +1,6 @@
 from utils import get_path
 from tacs import elements, constitutive, functions
+from mphys import MPhysVariables
 import numpy as np
 
 def get():
@@ -25,5 +26,6 @@ def get():
     return {
         #'element_callback': element_callback,
         'problem_setup': problem_setup,
+        'coupling_loads': [MPhysVariables.Structures.Loads.AERODYNAMIC],
         'mesh_file': get_path('crm.bdf', 'setup')
         }
diff --git a/crm/static.py b/crm/static.py
index b8a37066196590396a914f68886f651eb98198b0..d37c5579bdce5378ed839a5341a8a42ab113c012 100644
--- a/crm/static.py
+++ b/crm/static.py
@@ -1,12 +1,12 @@
 from utils import parse_args
-from setup import pygeo_setup, dart_setup, tacs_setup
+from setup import pygeo_setup, dart_setup, tacs_setup, meld_setup
 
 from dartflo.dart.api.mphys import DartBuilder
 from tacs.mphys import TacsBuilder
 from funtofem.mphys.mphys_meld import MeldBuilder
 from geometry import PyGeoBuilder
-from mphys.scenario_aerostructural import ScenarioAeroStructural
-from mphys import MultipointParallel
+from mphys.scenarios.aerostructural import ScenarioAeroStructural
+from mphys import MultipointParallel, MPhysVariables
 import openmdao.api as om
 
 class Static(MultipointParallel):
@@ -18,7 +18,7 @@ class Static(MultipointParallel):
             geom = PyGeoBuilder(pygeo_setup.get(i))
             dart = DartBuilder(dart_setup.get(args.k, args.v, i), scenario='aerostructural', task='optimization', saveAll=False, raiseError=False)
             tacs = TacsBuilder(**tacs_setup.get())
-            meld = MeldBuilder(dart, tacs, isym=1)
+            meld = MeldBuilder(dart, tacs, isym=1, **meld_setup.get())
             # create the scenarios
             csn = self.mphys_add_scenario(scn, ScenarioAeroStructural(aero_builder=dart, struct_builder=tacs, ldxfer_builder=meld, geometry_builder=geom, in_MultipointParallel=True), coupling_nonlinear_solver=om.NonlinearBlockGS(maxiter=10, iprint=2, use_aitken=True, rtol=1e-6, atol=1e-8, aitken_initial_factor=0.75), coupling_linear_solver=om.LinearBlockGS(maxiter=10, iprint=2, use_aitken=True, rtol=1e-6, atol=1e-8, aitken_initial_factor=0.75))
 
@@ -29,8 +29,8 @@ class Static(MultipointParallel):
         for i, ssys in enumerate(self.system_iter(recurse=False)):
             # set surface meshes
             ssys.geometry.setConstrainedSurface(ssys.aero_mesh.get_triangulated_surface())
-            ssys.geometry.addDisciplinePoints('aero', ssys.aero_mesh.get_val('x_aero0'))
-            ssys.geometry.addDisciplinePoints('struct', ssys.struct_mesh.masker.get_val('x_struct0_masked'))
+            ssys.geometry.addDisciplinePoints('aero', ssys.aero_mesh.get_val(MPhysVariables.Aerodynamics.Surface.Mesh.COORDINATES))
+            ssys.geometry.addDisciplinePoints('struct', ssys.struct_mesh.masker.get_val(f'{MPhysVariables.Structures.Mesh.COORDINATES}_masked'))
             # save sizes
             sizes = ssys.geometry.getSizes()
             self.geo_ntwist = sizes['twist']
diff --git a/m6/geometry.py b/m6/geometry.py
index e670214b06523a4ebd6352bf19f67f6ff47631fb..caf055557a6fe1d5b849dcc3cbedfb8eec9fbac2 100644
--- a/m6/geometry.py
+++ b/m6/geometry.py
@@ -4,7 +4,13 @@
 import numpy as np
 from mpi4py import MPI
 import openmdao.api as om
-from mphys.builder import Builder
+from mphys import Builder, MPhysVariables
+
+# Define short aliases for convenience
+XAERO_IN = MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_INPUT
+XAERO_OUT = MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_OUTPUT
+XSTRUCT_IN = MPhysVariables.Structures.Geometry.COORDINATES_INPUT
+XSTRUCT_OUT = MPhysVariables.Structures.Geometry.COORDINATES_OUTPUT
 
 # Geometry parametrization and constraints
 class PyGeoGeom(om.ExplicitComponent):
@@ -20,9 +26,12 @@ class PyGeoGeom(om.ExplicitComponent):
         self.omPtSetList = []
         # Setup geometry, constraints and their I/O
         # add IO for surface meshes, the DVGeo object must be configured later by calling addDisciplinePoints
-        for disc in self.config['disciplines']:
-            self.add_input(f'x_{disc}_in', distributed=True, shape_by_conn=True)
-            self.add_output(f'x_{disc}0', distributed=True, copy_shape=f'x_{disc}_in', tags=['mphys_coordinates'])
+        if 'aero' in self.config['disciplines']:
+            self.add_input(XAERO_IN, distributed=True, shape_by_conn=True, tags=['mphys_coordinates'])
+            self.add_output(XAERO_OUT, distributed=True, copy_shape=XAERO_IN, tags=['mphys_coordinates'])
+        if 'struct' in self.config['disciplines']:
+            self.add_input(XSTRUCT_IN, distributed=True, shape_by_conn=True, tags=['mphys_coordinates'])
+            self.add_output(XSTRUCT_OUT, distributed=True, copy_shape=XSTRUCT_IN, tags=['mphys_coordinates'])
         # add a reference axis
         if 'refAxis' in self.config:
             for name, prm in self.config['refAxis'].items():
@@ -88,8 +97,12 @@ class PyGeoGeom(om.ExplicitComponent):
     def addDisciplinePoints(self, discipline, points, **kwargs):
         """Add the points of the discipline to be embedded in the parametrization
         """
-        self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), f'x_{discipline}0', **kwargs)
-        self.omPtSetList.append(f'x_{discipline}0')
+        if discipline == 'aero':
+            self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), XAERO_OUT, **kwargs)
+            self.omPtSetList.append(XAERO_OUT)
+        elif discipline == 'struct':
+            self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), XSTRUCT_OUT, **kwargs)
+            self.omPtSetList.append(XSTRUCT_OUT)
 
     def setConstrainedSurface(self, surface):
         """Set (a triangulation of) the surface to be constrained and add the constraints to DVCon
@@ -291,6 +304,7 @@ class PyGeoBuilder(Builder):
         if self.__cfg['format'] == '.xyz':
             from pygeo import DVGeometry
             self.__dvGeo = DVGeometry(self.__cfg['geoFile'])
+            self.__dvGeo.writeTecplot('ffd_init.dat')
         else:
             from pygeo import DVGeometryVSP
             self.__dvGeo = DVGeometryVSP(self.__cfg['geoFile'], comm=self.comm, **self.__opts)
@@ -299,7 +313,7 @@ class PyGeoBuilder(Builder):
         self.__dvCon = DVConstraints()
         self.__dvCon.setDVGeo(self.__dvGeo)
 
-    def get_mesh_coordinate_subsystem(self):
+    def get_mesh_coordinate_subsystem(self, scenario_name=None):
         """Return OpenMDAO component that wraps pyGeo
         """
         return PyGeoGeom(dvGeo=self.__dvGeo, dvCon=self.__dvCon, config=self.__cfg)
diff --git a/m6/run.py b/m6/run.py
index 65ca13c55a17b8de7cccaec939be0a1aa52fe68d..f2b65575ceb4f431c694403f970a4d4f1f226775 100644
--- a/m6/run.py
+++ b/m6/run.py
@@ -4,8 +4,8 @@ from utils import change_workdir, parse_args, print_time
 from geometry import PyGeoBuilder
 from dartflo.dart.api.mphys import DartBuilder
 
-from mphys import Multipoint
-from mphys.scenario_aerodynamic import ScenarioAerodynamic
+from mphys import Multipoint, MPhysVariables
+from mphys.scenarios.aerodynamic import ScenarioAerodynamic
 import openmdao.api as om
 
 class Top(Multipoint):
@@ -28,8 +28,8 @@ class Top(Multipoint):
         # create the scenario
         self.mphys_add_scenario('cruise', ScenarioAerodynamic(aero_builder=dart))
         # connect
-        self.connect('mesh.x_aero0', 'geometry.x_aero_in')
-        self.connect('geometry.x_aero0', 'cruise.x_aero')
+        self.connect(f'mesh.{MPhysVariables.Aerodynamics.Surface.Mesh.COORDINATES}', f'geometry.{MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_INPUT}')
+        self.connect(f'geometry.{MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_OUTPUT}', f'cruise.{MPhysVariables.Aerodynamics.Surface.COORDINATES}')
 
     def configure(self):
         import numpy as np
@@ -38,7 +38,7 @@ class Top(Multipoint):
         self.connect('aoa', 'cruise.aoa')
         self.add_design_var('aoa', lower=0.0, upper=5.0, scaler=0.1, units='deg')
         # DV - shape
-        self.geometry.addDisciplinePoints('aero', self.mesh.get_val('x_aero0'))
+        self.geometry.addDisciplinePoints('aero', self.mesh.get_val(MPhysVariables.Aerodynamics.Surface.Mesh.COORDINATES))
         sizes = self.geometry.getSizes()
         self.dvs.add_output('twist', val=np.zeros(sizes['twist']))
         self.dvs.add_output('shape', val=np.zeros(sizes['shape']))
diff --git a/n12/geometry.py b/n12/geometry.py
index e670214b06523a4ebd6352bf19f67f6ff47631fb..caf055557a6fe1d5b849dcc3cbedfb8eec9fbac2 100644
--- a/n12/geometry.py
+++ b/n12/geometry.py
@@ -4,7 +4,13 @@
 import numpy as np
 from mpi4py import MPI
 import openmdao.api as om
-from mphys.builder import Builder
+from mphys import Builder, MPhysVariables
+
+# Define short aliases for convenience
+XAERO_IN = MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_INPUT
+XAERO_OUT = MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_OUTPUT
+XSTRUCT_IN = MPhysVariables.Structures.Geometry.COORDINATES_INPUT
+XSTRUCT_OUT = MPhysVariables.Structures.Geometry.COORDINATES_OUTPUT
 
 # Geometry parametrization and constraints
 class PyGeoGeom(om.ExplicitComponent):
@@ -20,9 +26,12 @@ class PyGeoGeom(om.ExplicitComponent):
         self.omPtSetList = []
         # Setup geometry, constraints and their I/O
         # add IO for surface meshes, the DVGeo object must be configured later by calling addDisciplinePoints
-        for disc in self.config['disciplines']:
-            self.add_input(f'x_{disc}_in', distributed=True, shape_by_conn=True)
-            self.add_output(f'x_{disc}0', distributed=True, copy_shape=f'x_{disc}_in', tags=['mphys_coordinates'])
+        if 'aero' in self.config['disciplines']:
+            self.add_input(XAERO_IN, distributed=True, shape_by_conn=True, tags=['mphys_coordinates'])
+            self.add_output(XAERO_OUT, distributed=True, copy_shape=XAERO_IN, tags=['mphys_coordinates'])
+        if 'struct' in self.config['disciplines']:
+            self.add_input(XSTRUCT_IN, distributed=True, shape_by_conn=True, tags=['mphys_coordinates'])
+            self.add_output(XSTRUCT_OUT, distributed=True, copy_shape=XSTRUCT_IN, tags=['mphys_coordinates'])
         # add a reference axis
         if 'refAxis' in self.config:
             for name, prm in self.config['refAxis'].items():
@@ -88,8 +97,12 @@ class PyGeoGeom(om.ExplicitComponent):
     def addDisciplinePoints(self, discipline, points, **kwargs):
         """Add the points of the discipline to be embedded in the parametrization
         """
-        self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), f'x_{discipline}0', **kwargs)
-        self.omPtSetList.append(f'x_{discipline}0')
+        if discipline == 'aero':
+            self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), XAERO_OUT, **kwargs)
+            self.omPtSetList.append(XAERO_OUT)
+        elif discipline == 'struct':
+            self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), XSTRUCT_OUT, **kwargs)
+            self.omPtSetList.append(XSTRUCT_OUT)
 
     def setConstrainedSurface(self, surface):
         """Set (a triangulation of) the surface to be constrained and add the constraints to DVCon
@@ -291,6 +304,7 @@ class PyGeoBuilder(Builder):
         if self.__cfg['format'] == '.xyz':
             from pygeo import DVGeometry
             self.__dvGeo = DVGeometry(self.__cfg['geoFile'])
+            self.__dvGeo.writeTecplot('ffd_init.dat')
         else:
             from pygeo import DVGeometryVSP
             self.__dvGeo = DVGeometryVSP(self.__cfg['geoFile'], comm=self.comm, **self.__opts)
@@ -299,7 +313,7 @@ class PyGeoBuilder(Builder):
         self.__dvCon = DVConstraints()
         self.__dvCon.setDVGeo(self.__dvGeo)
 
-    def get_mesh_coordinate_subsystem(self):
+    def get_mesh_coordinate_subsystem(self, scenario_name=None):
         """Return OpenMDAO component that wraps pyGeo
         """
         return PyGeoGeom(dvGeo=self.__dvGeo, dvCon=self.__dvCon, config=self.__cfg)
diff --git a/n12/run.py b/n12/run.py
index 30432895b01f9605c3bb51ac86ed133234d62e9d..da8ac236818293dc1808d114f4ccd7a12a7a1629 100644
--- a/n12/run.py
+++ b/n12/run.py
@@ -4,8 +4,8 @@ from utils import change_workdir, parse_args, print_time
 from geometry import PyGeoBuilder
 from dartflo.dart.api.mphys import DartBuilder
 
-from mphys import Multipoint
-from mphys.scenario_aerodynamic import ScenarioAerodynamic
+from mphys import Multipoint, MPhysVariables
+from mphys.scenarios.aerodynamic import ScenarioAerodynamic
 import openmdao.api as om
 
 class Top(Multipoint):
@@ -28,8 +28,8 @@ class Top(Multipoint):
         # create the scenario
         self.mphys_add_scenario('cruise', ScenarioAerodynamic(aero_builder=dart))
         # connect
-        self.connect('mesh.x_aero0', 'geometry.x_aero_in')
-        self.connect('geometry.x_aero0', 'cruise.x_aero')
+        self.connect(f'mesh.{MPhysVariables.Aerodynamics.Surface.Mesh.COORDINATES}', f'geometry.{MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_INPUT}')
+        self.connect(f'geometry.{MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_OUTPUT}', f'cruise.{MPhysVariables.Aerodynamics.Surface.COORDINATES}')
 
     def configure(self):
         import numpy as np
@@ -38,7 +38,7 @@ class Top(Multipoint):
         self.connect('aoa', 'cruise.aoa')
         self.add_design_var('aoa', lower=0.0, upper=5.0, scaler=0.1, units='deg')
         # DV - shape
-        self.geometry.addDisciplinePoints('aero', self.mesh.get_val('x_aero0'))
+        self.geometry.addDisciplinePoints('aero', self.mesh.get_val(MPhysVariables.Aerodynamics.Surface.Mesh.COORDINATES))
         sizes = self.geometry.getSizes()
         self.dvs.add_output('shape', val=np.zeros(sizes['shape']))
         self.connect('shape', 'geometry.shape')
diff --git a/rae_alu/geometry.py b/rae_alu/geometry.py
index a5905b7eac642ce2f5f2383ddbd1a2a0dfddf043..a0b75333035890b3a8f4e1ab2dec226405232328 100644
--- a/rae_alu/geometry.py
+++ b/rae_alu/geometry.py
@@ -5,7 +5,13 @@ import numpy as np
 from mpi4py import MPI
 import openmdao.api as om
 from omflut.builder import Builder as OMFlutBuilder
-from mphys.builder import Builder as MphysBuilder
+from mphys import Builder as MphysBuilder, MPhysVariables
+
+# Define short aliases for convenience
+XAERO_IN = MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_INPUT
+XAERO_OUT = MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_OUTPUT
+XSTRUCT_IN = MPhysVariables.Structures.Geometry.COORDINATES_INPUT
+XSTRUCT_OUT = MPhysVariables.Structures.Geometry.COORDINATES_OUTPUT
 
 # Geometry parametrization and constraints
 class PyGeoGeom(om.ExplicitComponent):
@@ -13,6 +19,7 @@ class PyGeoGeom(om.ExplicitComponent):
         self.options.declare('dvGeo', desc='geometry parametrization', recordable=False)
         self.options.declare('dvCon', desc='constraints', recordable=False)
         self.options.declare('config', desc='configuration', recordable=False)
+        self.options.declare('use_old_io', default=False, desc='use old in/out names', recordable=False)
 
     def setup(self):
         self.DVGeo = self.options['dvGeo']
@@ -21,9 +28,17 @@ class PyGeoGeom(om.ExplicitComponent):
         self.omPtSetList = []
         # Setup geometry, constraints and their I/O
         # add IO for surface meshes, the DVGeo object must be configured later by calling addDisciplinePoints
-        for disc in self.config['disciplines']:
-            self.add_input(f'x_{disc}_in', distributed=True, shape_by_conn=True)
-            self.add_output(f'x_{disc}0', distributed=True, copy_shape=f'x_{disc}_in', tags=['mphys_coordinates'])
+        if self.options['use_old_io']:
+            for disc in self.config['disciplines']:
+                self.add_input(f'x_{disc}_in', distributed=True, shape_by_conn=True)
+                self.add_output(f'x_{disc}0', distributed=True, copy_shape=f'x_{disc}_in', tags=['mphys_coordinates'])
+        else:
+            if 'aero' in self.config['disciplines']:
+                self.add_input(XAERO_IN, distributed=True, shape_by_conn=True, tags=['mphys_coordinates'])
+                self.add_output(XAERO_OUT, distributed=True, copy_shape=XAERO_IN, tags=['mphys_coordinates'])
+            if 'struct' in self.config['disciplines']:
+                self.add_input(XSTRUCT_IN, distributed=True, shape_by_conn=True, tags=['mphys_coordinates'])
+                self.add_output(XSTRUCT_OUT, distributed=True, copy_shape=XSTRUCT_IN, tags=['mphys_coordinates'])
         # add a reference axis
         if 'refAxis' in self.config:
             for name, prm in self.config['refAxis'].items():
@@ -96,8 +111,16 @@ class PyGeoGeom(om.ExplicitComponent):
     def addDisciplinePoints(self, discipline, points, **kwargs):
         """Add the points of the discipline to be embedded in the parametrization
         """
-        self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), f'x_{discipline}0', **kwargs)
-        self.omPtSetList.append(f'x_{discipline}0')
+        if self.options['use_old_io']:
+            self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), f'x_{discipline}0', **kwargs)
+            self.omPtSetList.append(f'x_{discipline}0')
+        else:
+            if discipline == 'aero':
+                self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), XAERO_OUT, **kwargs)
+                self.omPtSetList.append(XAERO_OUT)
+            elif discipline == 'struct':
+                self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), XSTRUCT_OUT, **kwargs)
+                self.omPtSetList.append(XSTRUCT_OUT)
 
     def setConstrainedSurface(self, surface):
         """Set (a triangulation of) the surface to be constrained and add the constraints to DVCon
@@ -352,4 +375,4 @@ class PyGeoBuilderOMF(OMFlutBuilder):
     def get_mesh(self):
         """Return OpenMDAO component that wraps pyGeo
         """
-        return PyGeoGeom(dvGeo=self.__dvGeo, dvCon=self.__dvCon, config=self.__cfg)
+        return PyGeoGeom(dvGeo=self.__dvGeo, dvCon=self.__dvCon, config=self.__cfg, use_old_io=True)
diff --git a/rae_alu/setup/meld_setup.py b/rae_alu/setup/meld_setup.py
new file mode 100644
index 0000000000000000000000000000000000000000..184a06ae7c029c1243c2f81ebfa764bf3d112077
--- /dev/null
+++ b/rae_alu/setup/meld_setup.py
@@ -0,0 +1,7 @@
+def get():
+    return {
+        'body_tags': [{
+            'aero': None,
+            'struct': ['SKIN', 'SPARS']
+            }]
+    }
diff --git a/rae_alu/setup/tacs_setup.py b/rae_alu/setup/tacs_setup.py
index f52ae0c73dcc8b8bdf64a8f78cfa0f26caca6c74..fb1bb69eb04b3ae491c0f7d8344f04aee047688b 100644
--- a/rae_alu/setup/tacs_setup.py
+++ b/rae_alu/setup/tacs_setup.py
@@ -1,5 +1,6 @@
 from utils import get_path
 from tacs import elements, constitutive, functions
+from mphys import MPhysVariables
 import numpy as np
 
 def get(pbl='static'):
@@ -82,6 +83,7 @@ def get(pbl='static'):
             'element_callback': element_callback,
             'problem_setup': problem_setup,
             'constraint_setup': constraint_setup,
+            'coupling_loads': [MPhysVariables.Structures.Loads.AERODYNAMIC],
             'mesh_file': get_path('wingbox-L1-Order2.bdf', 'setup')
             }
     elif pbl == 'modal':
diff --git a/rae_alu/static.py b/rae_alu/static.py
index 969b283f095b7e50027e7a0c2eae6408ac880452..11482fbf45884dba46270acf16862b3228b2d1a9 100644
--- a/rae_alu/static.py
+++ b/rae_alu/static.py
@@ -1,14 +1,12 @@
 from utils import parse_args
-from setup import pygeo_setup
-from setup import dart_setup
-from setup import tacs_setup
+from setup import pygeo_setup, dart_setup, tacs_setup, meld_setup
 
 from dartflo.dart.api.mphys import DartBuilder
 from tacs.mphys import TacsBuilder
 from funtofem.mphys.mphys_meld import MeldBuilder
 from geometry import PyGeoBuilder
-from mphys.scenario_aerostructural import ScenarioAeroStructural
-from mphys import MultipointParallel
+from mphys.scenarios.aerostructural import ScenarioAeroStructural
+from mphys import MultipointParallel, MPhysVariables
 import openmdao.api as om
 
 class Static(MultipointParallel):
@@ -20,7 +18,7 @@ class Static(MultipointParallel):
             geom = PyGeoBuilder(pygeo_setup.get(i))
             dart = DartBuilder(dart_setup.get(args.k, args.v, i), scenario='aerostructural', task='optimization', saveAll=False, raiseError=False)
             tacs = TacsBuilder(**tacs_setup.get())
-            meld = MeldBuilder(dart, tacs, isym=1)
+            meld = MeldBuilder(dart, tacs, isym=1, **meld_setup.get())
             # create the scenarios
             csn = self.mphys_add_scenario(scn, ScenarioAeroStructural(aero_builder=dart, struct_builder=tacs, ldxfer_builder=meld, geometry_builder=geom, in_MultipointParallel=True), coupling_nonlinear_solver=om.NonlinearBlockGS(maxiter=10, iprint=2, use_aitken=True, rtol=1e-6, atol=1e-8, aitken_initial_factor=0.75), coupling_linear_solver=om.LinearBlockGS(maxiter=10, iprint=2, use_aitken=True, rtol=1e-6, atol=1e-8, aitken_initial_factor=0.75))
 
@@ -31,8 +29,8 @@ class Static(MultipointParallel):
         for i, ssys in enumerate(self.system_iter(recurse=False)):
             # set surface meshes
             ssys.geometry.setConstrainedSurface(ssys.aero_mesh.get_triangulated_surface())
-            ssys.geometry.addDisciplinePoints('aero', ssys.aero_mesh.get_val('x_aero0'))
-            ssys.geometry.addDisciplinePoints('struct', ssys.struct_mesh.fea_mesh.get_val('x_struct0'))
+            ssys.geometry.addDisciplinePoints('aero', ssys.aero_mesh.get_val(MPhysVariables.Aerodynamics.Surface.Mesh.COORDINATES))
+            ssys.geometry.addDisciplinePoints('struct', ssys.struct_mesh.masker.get_val(f'{MPhysVariables.Structures.Mesh.COORDINATES}_masked'))
             # save DV
             if i == 0:
                 # geometrical design variables (size)
diff --git a/rae_cfrp/geometry.py b/rae_cfrp/geometry.py
index a5905b7eac642ce2f5f2383ddbd1a2a0dfddf043..a0b75333035890b3a8f4e1ab2dec226405232328 100644
--- a/rae_cfrp/geometry.py
+++ b/rae_cfrp/geometry.py
@@ -5,7 +5,13 @@ import numpy as np
 from mpi4py import MPI
 import openmdao.api as om
 from omflut.builder import Builder as OMFlutBuilder
-from mphys.builder import Builder as MphysBuilder
+from mphys import Builder as MphysBuilder, MPhysVariables
+
+# Define short aliases for convenience
+XAERO_IN = MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_INPUT
+XAERO_OUT = MPhysVariables.Aerodynamics.Surface.Geometry.COORDINATES_OUTPUT
+XSTRUCT_IN = MPhysVariables.Structures.Geometry.COORDINATES_INPUT
+XSTRUCT_OUT = MPhysVariables.Structures.Geometry.COORDINATES_OUTPUT
 
 # Geometry parametrization and constraints
 class PyGeoGeom(om.ExplicitComponent):
@@ -13,6 +19,7 @@ class PyGeoGeom(om.ExplicitComponent):
         self.options.declare('dvGeo', desc='geometry parametrization', recordable=False)
         self.options.declare('dvCon', desc='constraints', recordable=False)
         self.options.declare('config', desc='configuration', recordable=False)
+        self.options.declare('use_old_io', default=False, desc='use old in/out names', recordable=False)
 
     def setup(self):
         self.DVGeo = self.options['dvGeo']
@@ -21,9 +28,17 @@ class PyGeoGeom(om.ExplicitComponent):
         self.omPtSetList = []
         # Setup geometry, constraints and their I/O
         # add IO for surface meshes, the DVGeo object must be configured later by calling addDisciplinePoints
-        for disc in self.config['disciplines']:
-            self.add_input(f'x_{disc}_in', distributed=True, shape_by_conn=True)
-            self.add_output(f'x_{disc}0', distributed=True, copy_shape=f'x_{disc}_in', tags=['mphys_coordinates'])
+        if self.options['use_old_io']:
+            for disc in self.config['disciplines']:
+                self.add_input(f'x_{disc}_in', distributed=True, shape_by_conn=True)
+                self.add_output(f'x_{disc}0', distributed=True, copy_shape=f'x_{disc}_in', tags=['mphys_coordinates'])
+        else:
+            if 'aero' in self.config['disciplines']:
+                self.add_input(XAERO_IN, distributed=True, shape_by_conn=True, tags=['mphys_coordinates'])
+                self.add_output(XAERO_OUT, distributed=True, copy_shape=XAERO_IN, tags=['mphys_coordinates'])
+            if 'struct' in self.config['disciplines']:
+                self.add_input(XSTRUCT_IN, distributed=True, shape_by_conn=True, tags=['mphys_coordinates'])
+                self.add_output(XSTRUCT_OUT, distributed=True, copy_shape=XSTRUCT_IN, tags=['mphys_coordinates'])
         # add a reference axis
         if 'refAxis' in self.config:
             for name, prm in self.config['refAxis'].items():
@@ -96,8 +111,16 @@ class PyGeoGeom(om.ExplicitComponent):
     def addDisciplinePoints(self, discipline, points, **kwargs):
         """Add the points of the discipline to be embedded in the parametrization
         """
-        self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), f'x_{discipline}0', **kwargs)
-        self.omPtSetList.append(f'x_{discipline}0')
+        if self.options['use_old_io']:
+            self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), f'x_{discipline}0', **kwargs)
+            self.omPtSetList.append(f'x_{discipline}0')
+        else:
+            if discipline == 'aero':
+                self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), XAERO_OUT, **kwargs)
+                self.omPtSetList.append(XAERO_OUT)
+            elif discipline == 'struct':
+                self.DVGeo.addPointSet(points.reshape(len(points) // 3, 3), XSTRUCT_OUT, **kwargs)
+                self.omPtSetList.append(XSTRUCT_OUT)
 
     def setConstrainedSurface(self, surface):
         """Set (a triangulation of) the surface to be constrained and add the constraints to DVCon
@@ -352,4 +375,4 @@ class PyGeoBuilderOMF(OMFlutBuilder):
     def get_mesh(self):
         """Return OpenMDAO component that wraps pyGeo
         """
-        return PyGeoGeom(dvGeo=self.__dvGeo, dvCon=self.__dvCon, config=self.__cfg)
+        return PyGeoGeom(dvGeo=self.__dvGeo, dvCon=self.__dvCon, config=self.__cfg, use_old_io=True)
diff --git a/rae_cfrp/setup/meld_setup.py b/rae_cfrp/setup/meld_setup.py
new file mode 100644
index 0000000000000000000000000000000000000000..184a06ae7c029c1243c2f81ebfa764bf3d112077
--- /dev/null
+++ b/rae_cfrp/setup/meld_setup.py
@@ -0,0 +1,7 @@
+def get():
+    return {
+        'body_tags': [{
+            'aero': None,
+            'struct': ['SKIN', 'SPARS']
+            }]
+    }
diff --git a/rae_cfrp/setup/tacs_setup.py b/rae_cfrp/setup/tacs_setup.py
index 646bafd778ae31a971d65ca76302b5429d6bcf01..271fe2e0c582ee5c964b0581cbc9165ab8c118e7 100644
--- a/rae_cfrp/setup/tacs_setup.py
+++ b/rae_cfrp/setup/tacs_setup.py
@@ -1,5 +1,6 @@
 from utils import get_path
 from tacs import elements, constitutive, functions
+from mphys import MPhysVariables
 import numpy as np
 
 def get(pbl='static'):
@@ -105,6 +106,7 @@ def get(pbl='static'):
             'element_callback': element_callback,
             'problem_setup': problem_setup,
             'constraint_setup': constraint_setup,
+            'coupling_loads': [MPhysVariables.Structures.Loads.AERODYNAMIC],
             'mesh_file': get_path('wingbox-L1-Order2.bdf', 'setup')
             }
     elif pbl == 'modal':
diff --git a/rae_cfrp/static.py b/rae_cfrp/static.py
index c0c75fa459f48ad51d74188b8c49f5e1ff13dce4..9e3657c68bb12a0881af6b411ee1a58a497e8434 100644
--- a/rae_cfrp/static.py
+++ b/rae_cfrp/static.py
@@ -1,14 +1,12 @@
 from utils import parse_args
-from setup import pygeo_setup
-from setup import dart_setup
-from setup import tacs_setup
+from setup import pygeo_setup, dart_setup, tacs_setup, meld_setup
 
 from dartflo.dart.api.mphys import DartBuilder
 from tacs.mphys import TacsBuilder
 from funtofem.mphys.mphys_meld import MeldBuilder
 from geometry import PyGeoBuilder
-from mphys.scenario_aerostructural import ScenarioAeroStructural
-from mphys import MultipointParallel
+from mphys.scenarios.aerostructural import ScenarioAeroStructural
+from mphys import MultipointParallel, MPhysVariables
 import openmdao.api as om
 
 class Static(MultipointParallel):
@@ -20,7 +18,7 @@ class Static(MultipointParallel):
             geom = PyGeoBuilder(pygeo_setup.get(i))
             dart = DartBuilder(dart_setup.get(args.k, args.v, i), scenario='aerostructural', task='optimization', saveAll=False, raiseError=False)
             tacs = TacsBuilder(**tacs_setup.get())
-            meld = MeldBuilder(dart, tacs, isym=1)
+            meld = MeldBuilder(dart, tacs, isym=1, **meld_setup.get())
             # create the scenarios
             csn = self.mphys_add_scenario(scn, ScenarioAeroStructural(aero_builder=dart, struct_builder=tacs, ldxfer_builder=meld, geometry_builder=geom, in_MultipointParallel=True), coupling_nonlinear_solver=om.NonlinearBlockGS(maxiter=10, iprint=2, use_aitken=True, rtol=1e-6, atol=1e-8, aitken_initial_factor=0.75), coupling_linear_solver=om.LinearBlockGS(maxiter=10, iprint=2, use_aitken=True, rtol=1e-6, atol=1e-8, aitken_initial_factor=0.75))
 
@@ -32,8 +30,8 @@ class Static(MultipointParallel):
         for i, ssys in enumerate(self.system_iter(recurse=False)):
             # set surface meshes
             ssys.geometry.setConstrainedSurface(ssys.aero_mesh.get_triangulated_surface())
-            ssys.geometry.addDisciplinePoints('aero', ssys.aero_mesh.get_val('x_aero0'))
-            ssys.geometry.addDisciplinePoints('struct', ssys.struct_mesh.fea_mesh.get_val('x_struct0'))
+            ssys.geometry.addDisciplinePoints('aero', ssys.aero_mesh.get_val(MPhysVariables.Aerodynamics.Surface.Mesh.COORDINATES))
+            ssys.geometry.addDisciplinePoints('struct', ssys.struct_mesh.masker.get_val(f'{MPhysVariables.Structures.Mesh.COORDINATES}_masked'))
             # save DV
             if i == 0:
                 # geometrical design variables (size)