Merge branch 'adri' into 'master'

Version 1.1

See merge request !1

.gitlab-ci.yml

 default:
-    image: rboman/waves-py3:2020.3
+    image: rboman/waves-py3:2022.0
     before_script:
         # Instal pip
         - wget -q https://bootstrap.pypa.io/get-pip.py
@@ -10,6 +10,8 @@ default:
         - git checkout tags/v2.1.0
         - INCLUDE=${INCLUDE}:${PWD}/include
         - cd ..
+        # Install TACS
+        # ...
         # Install SDPM
         # TODO change branch, enable CoDiPack
         - wget -q https://gitlab.uliege.be/am-dept/sdpm/-/archive/adri/sdpm-adri.tar.bz2
@@ -22,15 +24,16 @@ default:
         - make install -j8
         - cd ../..
         # Install PyPK
-        - wget -q https://gitlab.uliege.be/am-dept/pypk/-/archive/master/pypk-master.tar.bz2
-        - tar xf pypk-master.tar.bz2
-        - rm pypk-master.tar.bz2
-        - mv pypk-master pypk
+        # TODO change branch
+        - wget -q https://gitlab.uliege.be/am-dept/pypk/-/archive/adri/pypk-adri.tar.bz2
+        - tar xf pypk-adri.tar.bz2
+        - rm pypk-adri.tar.bz2
+        - mv pypk-adri pypk
         - cd pypk
         - python3 -m pip install .
         - cd ..
         # Install OpenMDAO
-        - python3 -m pip install openmdao
+        - python3 -m pip install 'openmdao<3.35'
 
 .global_tag: &global_tag_def
     tags:

README.md

@@ -2,4 +2,4 @@
 OMFlut is a collection of tools to perform flutter analysis and flutter-constrained optimization calculations using [OpenMDAO](https://openmdao.org/). The package is written in Python and developed at the University of Liège.
 
 ## Documentation
-Detailed build and use instructions can be found in the [wiki](https://gitlab.uliege.be/am-dept/pypk/wikis/home).
+Detailed build and use instructions can be found in the [wiki](https://gitlab.uliege.be/am-dept/omflut/wikis/home).

examples/sdpm_nipk_agard.py

@@ -25,7 +25,7 @@ import numpy as np
 
 def get_path(fname):
     import os.path
-    return os.path.join(os.path.abspath(os.path.dirname(__file__)), 'data', os.path.basename(__file__)[:-3], fname)
+    return os.path.join(os.path.abspath(os.path.dirname(__file__)), 'data', fname)
 
 def get_exp():
     mach = [0.499, 0.678, 0.901, 0.960] # experimental Mach number
@@ -38,7 +38,7 @@ def get_cfg(irun):
     import os.path
     # Shared parameters
     c_ref = 0.559 # root chord
-    k_ref = [0.01, 0.1, 0.2, 0.3]#, 0.5, 0.7, 1.0] # reference reduced frequencies
+    k_ref = [0.01, 0.1, 0.2, 0.3] # reference reduced frequencies
     n_mod = 4 # number of mode
     # Experimental Mach number, density and mass ratio
     _, mach, rho, mu = get_exp()
@@ -51,6 +51,8 @@ def get_cfg(irun):
     # SDPM config
     wnc = 20
     a_cfg = {
+        # Verbosity
+        'Verb_lvl': 0,
         # Model (geometry or mesh)
         'File': get_path('agard445.geo'), # input Gmsh file
         'Pars': {'nC': wnc, 'nS': 10, 'pC': 1.1, 'pS': 1.0}, # parameters for Gmsh file (optional)
@@ -81,6 +83,7 @@ def get_cfg(irun):
     }
     # NIPK
     f_cfg = {
+        'name': 'agard', # case name
         'k_ref': np.array(k_ref), # reference reduced frequencies
         'l_ref': 0.5 * c_ref, # reference length (half root chord)
         'mach': mach[irun], # freesteam Mach number
@@ -103,7 +106,6 @@ class Agard(om.Group):
         self.options.declare('n', default=1, desc='number of runs')
 
     def setup(self):
-        self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*'])
         for i in range(self.options['n']):
             # Builders
             s_cfg, a_cfg, x_cfg, f_cfg = get_cfg(i)
@@ -114,15 +116,6 @@ class Agard(om.Group):
             # Flutter component
             self.add_subsystem(f'flutter_{i}', FlutterGroup(struct=agrd, xfer=rbfi, aero=sdpm, flutter=pypk))
 
-    def configure(self):
-        # DV
-        self.dvs.add_output('dummy', val=1, desc='dummy')
-        #self.add_design_var('dummy', lower=0.0, upper=1.0, scaler=1)
-        # CON
-        #self.add_constraint('flutter.ks_g', upper=-0.001, scaler=1e2)
-        # OBJ
-        #self.add_objective('flutter.m', scaler=1e-2)
-
 def print_time(selapsed):
     days, rem = divmod(selapsed, 24*60*60)
     hours, rem = divmod(rem, 60*60)

examples/tacs_sdpm_nipk_agard.py

+#!/usr/bin/env python3
+# -*- coding: utf8 -*-
+
+# Copyright 2024 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.
+
+from omflut import FlutterGroup
+from omflut.tacs_structure import TacsBuilder
+from omflut import RbfXferBuilder
+from sdpm.api.om import SdpmBuilder
+from pypk.api_om import PypkBuilder
+import openmdao.api as om
+import numpy as np
+
+def get_path(fname):
+    import os.path
+    return os.path.join(os.path.abspath(os.path.dirname(__file__)), 'data', fname)
+
+def get_exp():
+    mach = [0.499, 0.678, 0.901, 0.960] # experimental Mach number
+    rho = [0.4278, 0.2082, 0.0995, 0.0634] # experimental density
+    mu = [33.465, 68.753, 143.920, 225.820] # experimental mass ratio m/(rho*V)
+    mach = [0.960] # experimental Mach number
+    rho = [0.0634] # experimental density
+    mu = [225.820] # experimental mass ratio m/(rho*V)
+    n = len(mach)
+    return n, mach, rho, mu
+
+def get_cfg(irun):
+    # Shared parameters
+    c_ref = 0.559 # root chord
+    k_ref = [0.01, 0.1, 0.2, 0.3] # reference reduced frequencies
+    n_mod = 4 # number of modes
+    # Experimental Mach number, density and mass ratio
+    _, mach, rho, mu = get_exp()
+    # Structural config
+    s_cfg = tacs_cfg(n_mod)
+    # SDPM config
+    a_cfg = sdpm_cfg(c_ref, n_mod, mach[irun], k_ref)
+    # RBF
+    x_cfg = rbfi_cfg(n_mod)
+    # NIPK
+    f_cfg = pypk_cfg(c_ref, n_mod, mach[irun], rho[irun], mu[irun], k_ref)
+
+    return s_cfg, a_cfg, x_cfg, f_cfg
+
+def tacs_cfg(n_mod):
+    from tacs import constitutive, elements
+    # Material properties
+    rho = 381.9  # density, kg/m^3
+    e1 = 3.245e9 # Youngs' modulus in direction '1', Pa
+    e2 = 4.162e8 # Youngs' modulus in direction '2', Pa
+    nu12 = 0.31 # Poisson's ratio in plane '1-2'
+    g12 = 4.392e8 # shear modulus in plane '1-2', Pa
+    min_thickness = 0.0001
+    max_thickness = 0.05
+
+    # Element definition
+    def element_callback(dvNum, compID, compDescript, elemDescripts, specialDVs, **kwargs):
+        # Get the thickness of current element
+        thickness = float(compDescript.split('/')[1].split('=')[1])
+        # Set up orthotropic material and shell constitutive properties with 1 DV per element
+        prop = constitutive.MaterialProperties(rho=rho, E1=e1, E2=e2, nu12=nu12, G12=g12)
+        con = constitutive.IsoShellConstitutive(prop, t=thickness, tNum=dvNum, tlb=min_thickness, tub=max_thickness)
+        # Add the element
+        transform = elements.ShellRefAxisTransform(np.array([np.sqrt(2), np.sqrt(2), 0.0]))
+        return elements.Quad4Shell(transform, con)
+
+    # Various pyTACS options
+    options = {
+        'printtiming': True
+    }
+
+    return {
+        'mesh_file': get_path('agard445.bdf'),
+        'element_callback': element_callback,
+        'num_modes': n_mod,
+        'sigma': 1.,
+        'pytacs_options': options,
+        'write_solution': False
+        }
+
+def sdpm_cfg(c_ref, n_mod, mach, k_ref):
+    wnc = 20
+    return {
+        # Verbosity
+        'Verb_lvl': 0,
+        # Model (geometry or mesh)
+        'File': get_path('agard445.geo'), # input Gmsh file
+        'Pars': {'nC': wnc, 'nS': 10, 'pC': 1.1, 'pS': 1.0}, # parameters for Gmsh file (optional)
+        # Markers
+        'Wing': 'wing', # name of the marker defining the lifting body (wing)
+        'Te': 'te', # name of the marker defining the trailing edge of the lifting body
+        # Freestream
+        'Mach': mach, # Mach number (optional)
+        'AoA': 0.0, # angle of attack (optional)
+        'AoS': 0.0, # angle of side slip (optional)
+        # Unsteady motion
+        'Num_wake_div': wnc, # number of cells per chord length in the wake
+        'Frequencies': k_ref, # list of reduced frequencies
+        'Num_modes': n_mod, # number of modes
+        # Geometry
+        'Symmetry': True, # whether only half (symmetric) model is provided
+        's_ref': 0.353, # reference surface area
+        'c_ref': c_ref, # reference chord length
+        'x_ref': 0.0, # reference point for moment computation (x)
+        'y_ref': 0.0, # reference point for moment computation (y)
+        'z_ref': 0.0 # reference point for moment computation (z)
+    }
+
+def rbfi_cfg(n_mod):
+    return {
+    'n_dim': 2, # number of coordinates to use for interpolation
+    'n_modes': n_mod, # number of modes
+    'smooth': 1 # smoothing
+    }
+
+def pypk_cfg(c_ref, n_mod, mach, rho, mu, k_ref):
+    return {
+        'name': 'agard', # case name
+        'k_ref': np.array(k_ref), # reference reduced frequencies
+        'l_ref': 0.5 * c_ref, # reference length (half root chord)
+        'mach': mach, # freesteam Mach number
+        'rho_inf': rho, # freesteam density
+        'u_idx': np.linspace(0.2, 0.5, 31, endpoint=True), # velocity index range
+        'mu': mu, # mass ratio
+        'f_ref': 40.3511, # frequency of first torsional mode
+        'n_modes': n_mod, # number of modes
+        'rho_ks': 1e6, # KS aggregation parameter for KS
+        'method': 'nipk', # method type
+        'fluid': 'unmatched', # fluid type
+        'vrb': 1 # verbosity level
+    }
+
+def mass_builder(mesh_file, element_callback):
+    """Create TACS static problem to compute mass
+    """
+    from tacs.pytacs import pyTACS
+    from tacs import functions
+    from mpi4py import MPI
+    fea_assembler = pyTACS(mesh_file, comm=MPI.COMM_WORLD)
+    fea_assembler.initialize(element_callback)
+    pbl = fea_assembler.createStaticProblem('static')
+    pbl.addFunction('mass', functions.StructuralMass)
+    return pbl
+
+class AgardMass(om.ExplicitComponent):
+    """Wrap TACS static problem to compute wing mass
+    """
+    def initialize(self):
+        self.options.declare('pbl', desc='static problem', recordable=False)
+
+    def setup(self):
+        self.pbl = self.options['pbl']
+        self.add_input('dv_struct', shape_by_conn=True, desc='structural design variables')
+        self.add_output('m', val=0., desc='structural mass')
+        self.declare_partials(of=['m'], wrt=['dv_struct'], method='exact')
+
+    def compute(self, inputs, outputs):
+        self.pbl.setDesignVars(inputs['dv_struct'])
+        func = {}
+        self.pbl.evalFunctions(func, ['mass'])
+        outputs['m'] = func['static_mass']
+        self.pbl.writeSolution()
+
+    def compute_partials(self, inputs, partials):
+        func_sens = {}
+        self.pbl.evalFunctionsSens(func_sens, ['mass'])
+        partials['m', 'dv_struct'] = func_sens['static_mass']['struct']
+
+class Agard(om.Group):
+    """Aerostructural flutter model of the AGARD 445.6 wing
+    """
+    def initialize(self):
+        self.options.declare('n', default=1, desc='number of runs')
+
+    def setup(self):
+        self.n_run = self.options['n']
+        self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*'])
+        for i in range(self.n_run):
+            # Builders
+            s_cfg, a_cfg, x_cfg, f_cfg = get_cfg(i)
+            tacs = TacsBuilder(**s_cfg)
+            sdpm = SdpmBuilder(a_cfg)
+            rbfi = RbfXferBuilder(tacs, sdpm, x_cfg)
+            pypk = PypkBuilder(f_cfg)
+            # Flutter component
+            self.add_subsystem(f'flutter_{i}', FlutterGroup(struct=tacs, xfer=rbfi, aero=sdpm, flutter=pypk))
+        # Mass component
+        pbl = mass_builder(s_cfg['mesh_file'], s_cfg['element_callback'])
+        self.init_dvs = pbl.getDesignVars()
+        self.add_subsystem('mass', AgardMass(pbl=pbl), promotes=['dv_struct'])
+
+    def configure(self):
+        # DV
+        self.dvs.add_output('dv_struct', val=self.init_dvs, desc='structural DV')
+        self.add_design_var('dv_struct', lower=0.0001, upper=0.05, scaler=1e3)
+        for i in range(self.n_run):
+            self.connect('dv_struct', f'flutter_{i}.struct.dv_struct')
+        # CON
+        for i in range(self.n_run):
+            self.add_constraint(f'flutter_{i}.g_agg', upper=-0.001, scaler=1e3)
+        # OBJ
+        self.add_objective('mass.m', scaler=1)
+
+def print_time(selapsed):
+    days, rem = divmod(selapsed, 24*60*60)
+    hours, rem = divmod(rem, 60*60)
+    minutes, seconds = divmod(rem, 60)
+    print('Wall-clock time: {:0>2}-{:0>2}:{:0>2}:{:0>2}'.format(int(days),int(hours),int(minutes),int(seconds)))
+
+if __name__ == "__main__":
+    import time
+    # Set up the problem
+    n_runs, mach, _, _ = get_exp()
+    prob = om.Problem(model=Agard(n=n_runs), driver=om.ScipyOptimizeDriver(optimizer='SLSQP', tol=1e-4, disp=True, maxiter=20))
+    rec_inc = []
+    for i in range(n_runs):
+        rec_inc.extend([f'flutter_{i}.f', f'flutter_{i}.g', f'flutter_{i}.g_agg', f'flutter_{i}.flutter.f_speed', 'mass.m'])
+    prob.model.recording_options['includes'] = rec_inc
+    prob.model.add_recorder(om.SqliteRecorder('cases.sql'))
+    prob.setup()
+    om.n2(prob, show_browser=False, outfile='n2.html')
+
+    # Run
+    tic = time.perf_counter()
+    #prob.run_model()
+    #prob.check_partials(excludes=['*mass', '*aero', '*struct', '*xfer'], compact_print=False, form='central')
+    #prob.check_totals(of=['flutter_0.g_agg'], wrt=['dv_struct'], compact_print=False, form='central', show_progress=True, driver_scaling=False)
+    prob.run_driver()
+    toc = time.perf_counter()
+    print_time(toc-tic)
+
+    # Print the history
+    cases = om.CaseReader('cases.sql').get_cases()
+    hdr = '{:>10s}'.format('it')
+    for irun in range(n_runs):
+        hdr += ', {:>10s}, {:>10s}, {:>10s}'.format(f'g_agg_{irun}', f'u_f_{irun}', 'm')
+    print(hdr)
+    uf = np.zeros(n_runs)
+    for i, case in enumerate(cases):
+        val = '{:10d}'.format(i)
+        for irun in range(n_runs):
+            g_agg = case.get_val(f'flutter_{irun}.g_agg')[0]
+            uf[irun] = case.get_val(f'flutter_{irun}.flutter.f_speed')[0]
+            m = case.get_val('mass.m')[0]
+            val += ', {:10.4f}, {:10.4f}, {:10.4f}'.format(g_agg, uf[irun], m)
+        print(val)
+    # Save results
+    np.savetxt('mach_uf.dat', np.column_stack((mach, uf)), delimiter=', ', fmt='%15.2f', header='{:>13s}, {:>15s}'.format('Mach number', 'Flutter speed'))

omflut/__init__.py

@@ -15,7 +15,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-__version__ = '1.0.0'
+__version__ = '1.1.0'
 
 from .builder import Builder
 from .coupling import FlutterGroup

omflut/builder.py

@@ -27,7 +27,7 @@ class Builder():
         """
         raise NotImplementedError()
 
-    def get_solver(self):
+    def get_solver(self, scenario_name=''):
         """Return OpenMDAO component containing the solver
         """
         raise NotImplementedError()

omflut/coupling.py

@@ -24,11 +24,19 @@ class FlutterGroup(om.Group):
         self.options.declare('xfer', desc='interpolator', recordable=False)
         self.options.declare('aero', desc='aerodynamic solver', recordable=False)
         self.options.declare('flutter', desc='flutter solver', recordable=False)
+        self.options.declare('geo', desc='optional geometry parametrization', default=None, recordable=False)
 
     def setup(self):
-        self.add_subsystem('mesh_struct', self.options['struct'].get_mesh(), promotes=['x_struct0'])
-        self.add_subsystem('mesh_aero', self.options['aero'].get_mesh(), promotes=['x_aero0'])
-        self.add_subsystem('struct', self.options['struct'].get_solver(), promotes=['x_struct0', 'q_struct', 'M', 'K'])
-        self.add_subsystem('xfer', self.options['xfer'].get_solver(), promotes=['x_struct0', 'q_struct', 'x_aero0', 'q_aero'])
-        self.add_subsystem('aero', self.options['aero'].get_solver(), promotes=['x_aero0', 'q_aero', 'Q_re', 'Q_im'])
-        self.add_subsystem('flutter', self.options['flutter'].get_solver(), promotes=['M', 'K', 'Q_re', 'Q_im', 'f', 'g', 'g_agg'])
+        if self.options['geo'] is None:
+            self.add_subsystem('mesh_struct', self.options['struct'].get_mesh(), promotes=['x_struct0'])
+            self.add_subsystem('mesh_aero', self.options['aero'].get_mesh(), promotes=['x_aero0'])
+        else:
+            self.add_subsystem('mesh_struct', self.options['struct'].get_mesh())
+            self.add_subsystem('mesh_aero', self.options['aero'].get_mesh())
+            self.add_subsystem('geometry', self.options['geo'].get_mesh(), promotes=['x_struct0', 'x_aero0'])
+            self.connect('mesh_struct.x_struct0', 'geometry.x_struct_in')
+            self.connect('mesh_aero.x_aero0', 'geometry.x_aero_in')
+        self.add_subsystem('struct', self.options['struct'].get_solver(self.name), promotes=['x_struct0', 'q_struct', 'M', 'K'])
+        self.add_subsystem('xfer', self.options['xfer'].get_solver(self.name), promotes=['x_struct0', 'q_struct', 'x_aero0', 'q_aero'])
+        self.add_subsystem('aero', self.options['aero'].get_solver(self.name), promotes=['x_aero0', 'q_aero', 'Q_re', 'Q_im'])
+        self.add_subsystem('flutter', self.options['flutter'].get_solver(self.name), promotes=['M', 'K', 'Q_re', 'Q_im', 'f', 'g', 'g_agg'])

omflut/fake_structure.py

@@ -69,7 +69,7 @@ class FakeStructBuilder(Builder):
         """
         return FakeStructMesh(xs=self.__x_struct)
 
-    def get_solver(self):
+    def get_solver(self, scenario_name=''):
         """Return OpenMDAO component containing the solver
         """
         return FakeStructSolver(xs=self.__x_struct, qm=self.__q_modes, m=self.__m_mat, k=self.__k_mat)

omflut/rbf_xfer.py

@@ -43,7 +43,7 @@ class RbfXferBuilder(Builder):
         n_a = aero_builder.get_number_of_nodes() # number of aerodynamic nodes
         self.xfer = RbfInterp(cfg, n_s, n_a)
 
-    def get_solver(self):
+    def get_solver(self, scenario_name=''):
         """Return OpenMDAO component containing the solver
         """
         return RbfXferSolver(sol=self.xfer)

omflut/tacs_structure.py

+from omflut import Builder
+import openmdao.api as om
+import numpy as np
+
+class TacsMesh(om.IndepVarComp):
+    """Initial structural mesh
+    """
+    def initialize(self):
+        self.options.declare('fea_assembler', desc='pyTACS interface', recordable=False)
+        self.options.declare('x_mask', desc='array to recover true degrees of freedom (node coordinates)')
+
+    def setup(self):
+        asb = self.options['fea_assembler']
+        mask = self.options['x_mask']
+        # Store all the node coordinates and add as output
+        self.add_output('x_struct0', val=asb.getOrigNodes()[mask], desc='initial structural node coordinates')
+
+class TacsSolver(om.ExplicitComponent):
+    """OpenMDAO structural wrapper
+    
+    Attributes
+    ----------
+    n_modes : int
+        number of eigenvalue frequencies / eigenvector modes
+    abs : tacs.pyTACS object
+        pyTACS interface
+    pbl : tacs.problem.modal object
+        Modal problem definition
+    x_mask : np.array
+        Array to recover true degrees of freedom (node coordinates)
+    q_mask : np.array
+        Array to recover true degrees of freedom (mode shapes)
+    """
+    def initialize(self):
+        self.options.declare('fea_assembler', desc='pyTACS interface', recordable=False)
+        self.options.declare('fea_problem', desc='modal problem', recordable=False)
+        self.options.declare('write_solution', desc='flag to write solution')
+        self.options.declare('x_mask', desc='array to recover true degrees of freedom (node coordinates)')
+        self.options.declare('q_mask', desc='array to recover true degrees of freedom (mode shapes)')
+        self.options.declare('n_nodes', desc='number of retained nodes')
+        self.options.declare('retained_modes', desc='list of retained modes')
+
+    def setup(self):
+        self.asb = self.options['fea_assembler']
+        self.pbl = self.options['fea_problem']
+        self.x_mask = self.options['x_mask']
+        self.q_mask = self.options['q_mask']
+        self.retained_modes = self.options['retained_modes']
+        self.n_modes = len(self.retained_modes)
+        self.add_input('dv_struct', shape_by_conn=True, desc='structural design variables')
+        self.add_input('x_struct0', shape_by_conn=True, desc='structural coordinates')
+        self.add_output('q_struct', val=np.zeros((3 * self.options['n_nodes'], self.n_modes)), desc='modal displacements')
+        self.add_output('M', val=np.identity(self.n_modes), desc='mass matrix')
+        self.add_output('K', val=np.zeros((self.n_modes, self.n_modes)), desc='stiffness matrix')
+        # Partials
+        # self.declare_partials(of=['q_struct', 'M', 'K'], wrt=['x_struct0'], method='exact')
+        self.declare_partials(of=['K'], wrt=['dv_struct'], method='exact')
+
+    def compute(self, inputs, outputs):
+        # Update nodes and design variables
+        x_s = self.asb.getOrigNodes()
+        x_s[self.x_mask] = inputs['x_struct0']
+        self.pbl.setNodes(x_s)
+        self.pbl.setDesignVars(inputs['dv_struct'])
+        # Perform modal analysis
+        self.pbl.solve()
+        if self.options['write_solution']:
+            self.pbl.writeSolution()
+        # Mask and set modes
+        for i in range(self.n_modes):
+            _, q_s = self.pbl.getVariables(self.retained_modes[i])
+            outputs['q_struct'][:, i] = self._extract_data(q_s[self.q_mask])
+        # Get eigenvalues and set modal stiffness matrix
+        funcs = {}
+        self.pbl.evalFunctions(funcs)
+        outputs['K'] = np.diag(list(funcs.values()))[np.ix_(self.retained_modes, self.retained_modes)]
+
+    def compute_partials(self, inputs, partials):
+        # Approximate derivative of stiffness matrix diagonal entries by derivative of eigenvalues
+        func_sens = {}
+        self.pbl.evalFunctionsSens(func_sens, evalVars='dv')
+        for i in range(self.n_modes):
+            d_lambda = func_sens[f'{self.pbl.name}_eigsm.{self.retained_modes[i]}']
+            partials['K', 'dv_struct'][np.ravel_multi_index(([i], [i]), (self.n_modes, self.n_modes)), :] = d_lambda['struct']
+
+    def _extract_data(self, q):
+        """Extract modal displacements along z and rotations about x and y
+        """
+        # Constants
+        vars_node = 6 # number of variables per node
+        v_idx = [2, 3, 4] # indices of dz, rx and ry
+        n_var = len(v_idx) # number of variables to extract
+        n_nod = len(q) // vars_node # number of nodes
+        # Extract
+        q_red = np.zeros(n_nod * n_var)
+        for i in range(n_nod):
+            for j in range(n_var):
+                q_red[i * n_var + j] = q[i * vars_node + v_idx[j]]
+        return q_red
+
+class TacsBuilder(Builder):
+    """TACS builder for OpenMDAO
+
+    Attributes
+    ----------
+    fea_assembler : tacs.pyTACS object
+        pyTACS interface
+    fea_problem : tacs.problem.modal object
+        Modal problem definition
+    mask : np.array
+        Array to recover true degrees of freedom
+    write_solution : bool
+        Flag indicating whether to write solution
+    """
+    def __init__(self, mesh_file, num_modes, sigma=1., element_callback=None, exclude_modes=None, components_tag=None, pytacs_options=None, write_solution=True):
+        """Instantiate and initialize TACS components
+
+        Parameters
+        ----------
+        mesh_file : str or pyNastran.bdf.bdf.BDF
+            The BDF file or a pyNastran BDF object to load
+        num_modes : int
+            Number of modes
+        sigma : float
+            Guess for the lowest eigenvalue (default: 1.0)
+        element_callback : collections.abc.Callable, optional
+            User-defined callback function for setting up TACS elements and element DVs (default: None)
+        exclude_modes : list[int]
+            List of indices of modes to exclude from analysis (default: None)
+        components_tag : list[str]
+            List of tag names for which nodes will be included (default: None)
+        pytacs_options : dict, optional
+            Options dictionary passed to pyTACS assembler (default: None)
+        write_solution : bool, optional
+            Flag to determine whether to write out TACS solutions to f5 file each design iteration (default: True)
+        """
+        from tacs.pytacs import pyTACS
+        from mpi4py import MPI
+        # Initialize assembler
+        self.fea_assembler = pyTACS(mesh_file, options=pytacs_options, comm=MPI.COMM_WORLD)
+        self.fea_assembler.initialize(element_callback)
+        # Set up the problem
+        self.fea_problem = self.fea_assembler.createModalProblem('modal', sigma, num_modes)
+        self.retained_modes = list(range(num_modes))
+        if exclude_modes is not None: self.retained_modes = list(np.delete(self.retained_modes, exclude_modes))
+        # Select nodes from components
+        nodes_id = self._select_nodes(components_tag)
+        self.n_nodes = len(nodes_id)
+        # Create masks
+        self.x_mask = self._create_mask(nodes_id, 3)
+        self.q_mask = self._create_mask(nodes_id, 6)
+        # Save other parameters
+        self.write_solution = write_solution
+
+    def get_mesh(self):
+        """Return OpenMDAO component to get the initial mesh coordinates
+        """
+        return TacsMesh(fea_assembler=self.fea_assembler, x_mask=self.x_mask)
+
+    def get_solver(self, scenario_name=''):
+        """Return OpenMDAO component containing the solver
+        """
+        return TacsSolver(fea_assembler=self.fea_assembler, fea_problem=self.fea_problem, write_solution=self.write_solution, x_mask=self.x_mask, q_mask=self.q_mask, n_nodes=self.n_nodes, retained_modes=self.retained_modes)
+
+    def get_number_of_nodes(self):
+        """Return the number of (true) nodes
+        """
+        return self.n_nodes
+
+    def get_number_of_dv(self):
+        """Return the number of degrees of freedom
+        """
+        return self.fea_assembler.getTotalNumDesignVars()
+
+    def _select_nodes(self, components):
+        """Select nodes from given components
+        
+        Parameters
+        ----------
+        components : list[str]
+            List of tag names for which nodes will be included
+        """
+        if components is None:
+            nodes_id = np.arange(self.fea_assembler.getNumOwnedNodes())
+        else:
+            nodes_id = self.fea_assembler.getLocalNodeIDsForComps(self.fea_assembler.selectCompIDs(components))
+        mults_id = self.fea_assembler.getLocalMultiplierNodeIDs()
+        return list(set(nodes_id) - set(mults_id))
+
+    def _create_mask(self, nodes_id, vars_node):
+        """Create a mask to recover array indices associated to retained true degrees of freedom
+        
+        Parameters
+        ----------
+        nodes_id : list[int]
+            index of retained nodes
+        vars_node : int
+            Number of variables per node
+        """
+        mask = np.full((self.fea_assembler.getNumOwnedNodes(), vars_node), False)
+        mask[nodes_id, :] = True
+        return mask.flatten()