# reconstructSolution

Recover full-model transient solution from reduced-order model (ROM)

## Syntax

## Description

recovers the full transient thermal solution from the reduced-order model
`thermalresults`

= reconstructSolution(`Rtherm`

,`u_therm`

,`tlist`

)`Rtherm`

, temperature in modal coordinates
`u_therm`

, and the time-steps `tlist`

that you used to
solve the reduced model.

## Examples

### Reconstruct Structural Solution from ROM Results

Knowing the solution in terms of the interface degrees of freedom (DoFs) and modal DoFs, reconstruct the solution for the full structural transient analysis.

**Define Parameters for Structural Analysis**

Create a square cross-section beam geometry.

gm = multicuboid(0.05,0.003,0.003);

Plot the geometry, displaying face and edge labels.

```
figure
pdegplot(gm,FaceLabels="on",FaceAlpha=0.5)
view([71 4])
```

```
figure
pdegplot(gm,EdgeLabels="on",FaceAlpha=0.5)
view([71 4])
```

Add a vertex at the center of face 3.

centerVertex = addVertex(gm,Coordinates=[0.025 0.0 0.0015]);

Create an `femodel`

object for transient structural analysis and include the geometry in the model.

model = femodel(AnalysisType="structuralTransient", ... Geometry=gm);

Specify Young's modulus, Poisson's ratio, and the mass density of the material.

model.MaterialProperties = ... materialProperties(YoungsModulus=210E9, ... PoissonsRatio=0.3, ... MassDensity=7800);

Fix one end of the beam.

`model.EdgeBC([2 8 11 12]) = edgeBC(Constraint="fixed");`

Generate a mesh.

model = generateMesh(model);

Apply a sinusoidal concentrated force in the *z*-direction on the new vertex. First, define a sinusoidal load function, `sinusoidalLoad`

, to model a harmonic load. This function accepts the load magnitude (amplitude), `location`

and `state`

structure arrays, frequency, and phase. Because the function depends on time, it must return a matrix of `NaN`

of the correct size when `state.time`

is `NaN`

. Solvers check whether a problem is nonlinear or time-dependent by passing `NaN`

state values and looking for returned `NaN`

values.

function Tn = sinusoidalLoad(load,location,state,Frequency,Phase) if isnan(state.time) normal = [location.nx location.ny]; if isfield(location,"nz") normal = [normal location.nz]; end Tn = NaN*normal; return end if isa(load,"function_handle") load = load(location,state); else load = load(:); end % Transient model excited with harmonic load Tn = load.*sin(Frequency.*state.time + Phase); end

Now apply the force on the new vertex.

```
Force = [0 0 10];
Frequency = 6000;
Phase = 0;
forcePulse = @(location,state) ...
sinusoidalLoad(Force,location,state,Frequency,Phase);
model.VertexLoad(centerVertex) = vertexLoad(Force=forcePulse);
```

Specify zero initial conditions.

model.CellIC = cellIC(Velocity=[0 0 0],Displacement=[0 0 0]);

**Reduce Model**

Specify the fixed and loaded boundaries as structural superelement interfaces by creating a `romInterface`

object for each superelement interface. The reduced-order model technique retains the DoFs on the superelement interfaces while condensing all other DoFs to a set of modal DoFs. For better performance, use the set of edges bounding face 5 instead of using the entire face.

romObj1 = romInterface(Edge=[2 8 11 12]); romObj2 = romInterface(Vertex=centerVertex);

Assign a vector of interface objects to the `ROMInterfaces`

property of the model.

model.ROMInterfaces = [romObj1,romObj2];

Reduce the structure, retaining all fixed interface modes up to `5e5`

.

rom = reduce(model,FrequencyRange=[-0.1,5e5]);

**Simulate Transient Dynamics Using ROM**

Next, use the reduced-order model to simulate the transient dynamics. Use the `ode15s`

function directly to integrate the reduced system of ordinary differential equations. Take the loaded and modal DoFs for time-integration, and leave the fixed DoFs aside because the solution remains zero for those DoFs.

Working with the reduced model requires indexing into the reduced system matrices `rom.K`

and `rom.M`

. The arrangement of DoFs in reduced system is such that the physical DoFs corresponding to retained interfaces appear first followed by the generalized model DoFs. DoFs in a structural problem correspond to translational displacements. If the number of mesh points in a model is `Nn`

, then the software assigns the IDs to the DoFs as follows: the first `1`

to `Nn`

are *x*-displacements, `Nn+1`

to `2*Nn`

are *y*-displacements, and `2Nn+1`

to `3*Nn`

are *z*-displacements. Only the subset of these `3*Nn`

DoFs corresponding to `ROMInterfaces`

is retained in the reduced model. The reduced model object `rom`

contains these IDs for the retained DoFs in `rom.RetainedDoF`

.

Create a function that returns DoF IDs given node IDs and the number of nodes.

getDoF = @(x,numNodes) [x(:); x(:) + numNodes; x(:) + 2*numNodes];

Find the node at the loaded vertex.

`loadedNode = findNodes(rom.Mesh,"region",Vertex=centerVertex);`

Find the DoF of the loaded nodes using the helper function `getDoF`

.

numNodes = size(rom.Mesh.Nodes,2); loadDoFs = getDoF(loadedNode,numNodes);

Knowing the DoF IDs for the given node IDs, use `rom.RetainedDoF`

and the `intersect`

function to find the required indices corresponding to those DoF in the reduced matrices.

[~,loadNodeROMIds] = intersect(rom.RetainedDoF,loadDoFs);

In the reduced matrices `rom.K`

and `rom.M`

, generalized modal DoFs appear after the retained DoFs. Find the indices of modal DoFs in `rom`

matrices.

modelDoFIDs = ((numel(rom.RetainedDoF) + 1):size(rom.K,1))';

Find the indices for the ODE DoFs in reduced matrices. Because fixed-end DoFs are not a part of the ODE system, these indices are as follows.

odeDoFs = [loadNodeROMIds;modelDoFIDs];

Find the relevant components of `rom.K`

and `rom.M`

for time integration.

Kconstrained = rom.K(odeDoFs,odeDoFs); Mconstrained = rom.M(odeDoFs,odeDoFs); numODE = numel(odeDoFs);

Now you have a second-order system of ODEs. To use `ode15s`

, you must convert this system into a system of first-order ODEs by applying linearization. This type of a first-order system is twice the size of the second-order system.

Mode = [eye(numODE,numODE), zeros(numODE,numODE); ... zeros(numODE,numODE), Mconstrained]; Kode = [zeros(numODE,numODE), -eye(numODE,numODE); ... Kconstrained, zeros(numODE,numODE)]; Fode = zeros(2*numODE,1);

The specified concentrated force load in the full system is along the *z*-direction, which is the third DoF in the ODE system. Accounting for the linearization, obtain the first-order system to get the loaded ODE DoF.

loadODEDoF = numODE + 3;

Specify the mass matrix and the Jacobian for the ODE solver.

odeoptions = odeset; odeoptions = odeset(odeoptions,"Jacobian",-Kode); odeoptions = odeset(odeoptions,"Mass",Mode);

Specify zero initial conditions.

u0 = zeros(2*numODE,1);

Solve the reduced system by using ode15s and the helper function.

function f = CMSODEf(t,u,Kode,Fode,centerVertex) Fode(centerVertex) = 10*sin(6000*t); f = -Kode*u +Fode; end tlist = 0:0.00005:3E-3; sol = ode15s(@(t,y) CMSODEf(t,y,Kode,Fode,loadODEDoF), ... tlist,u0,odeoptions);

Compute the values of the ODE variable and the time derivatives.

[displ,vel] = deval(sol,tlist);

**Reconstruct Solution for Full Model**

Knowing the solution in terms of the interface DoFs and modal DoFs, you can reconstruct the solution for the full model. The `reconstructSolution`

function requires the displacement, velocity, and acceleration at all DoFs in `rom`

. Create the complete solution vector, including the zero values at the fixed DoFs.

u = zeros(size(rom.K,1),numel(tlist)); ut = zeros(size(rom.K,1),numel(tlist)); utt = zeros(size(rom.K,1),numel(tlist)); u(odeDoFs,:) = displ(1:numODE,:); ut(odeDoFs,:) = vel(1:numODE,:); utt(odeDoFs,:) = vel(numODE+1:2*numODE,:);

Create a transient results object using this solution.

RTrom = reconstructSolution(rom,u,ut,utt,tlist);

Compute the displacement in the interior at the center of the beam using the reconstructed solution.

```
coordCenter = [0;0;0];
iDispRTrom = interpolateDisplacement(RTrom,coordCenter);
figure
plot(tlist,iDispRTrom.uz)
title("Z-Displacement at Geometric Center")
```

### Reconstruct Thermal Solution from ROM Results

Reconstruct the solution for a full thermal transient analysis from the reduced-order model.

Create an `femodel`

object for transient thermal analysis, and include a unit square geometry in the model.

model = femodel(AnalysisType="thermalTransient", ... Geometry=@squareg);

Plot the geometry, displaying edge labels.

```
pdegplot(model,EdgeLabels="on")
xlim([-1.1 1.1])
ylim([-1.1 1.1])
```

Specify the thermal conductivity, mass density, and specific heat of the material.

model.MaterialProperties = ... materialProperties(ThermalConductivity=400, ... MassDensity=1300, ... SpecificHeat=600);

Set the temperature on the right edge to `100`

.

model.EdgeBC(2) = edgeBC(Temperature=100);

Set an initial value of 5`0`

for the temperature.

model.FaceIC = faceIC(Temperature=50);

Generate a mesh.

model = generateMesh(model);

Solve the model for three different values of heat source, and collect snapshots.

tlist = 0:10:600; snapShotIDs = [1:10 59 60 61]; Tmatrix = []; heatVariation = [10000 15000 20000]; for q = heatVariation model.FaceLoad = faceLoad(Heat=q); results = solve(model,tlist); Tmatrix = [Tmatrix,results.Temperature(:,snapShotIDs)]; end

Switch the thermal model analysis type to modal.

`model.AnalysisType = "thermalModal";`

Compute the POD modes.

RModal = solve(model,Snapshots=Tmatrix);

Reduce the thermal model.

Rtherm = reduce(model,ModalResults=RModal)

Rtherm = ReducedThermalModel with properties: K: [6x6 double] M: [6x6 double] F: [6x1 double] InitialConditions: [6x1 double] Mesh: [1x1 FEMesh] ModeShapes: [1529x5 double] SnapshotsAverage: [1529x1 double]

Next, use the reduced-order model to simulate the transient dynamics. Use the `ode15s`

function directly to integrate the reduced system ODE. Specify the mass matrix and the Jacobian for the ODE solver.

```
odeoptions = odeset;
odeoptions = odeset(odeoptions,Mass=Rtherm.M);
odeoptions = odeset(odeoptions,JConstant="on");
f = @(t,u) -Rtherm.K*u + Rtherm.F;
df = -Rtherm.K;
odeoptions = odeset(odeoptions,Jacobian=df);
```

Solve the reduced system by using `ode15s`

.

sol = ode15s(f,tlist,Rtherm.InitialConditions,odeoptions);

Compute the values of the ODE variable.

u = deval(sol,tlist);

Reconstruct the solution for the full model.

R = reconstructSolution(Rtherm,u,tlist);

Plot the temperature distribution at the last time step.

pdeplot(R.Mesh,XYData=R.Temperature(:,end))

## Input Arguments

`Rcb`

— Structural results obtained using Craig-Bampton order reduction method

`ReducedStructuralModel`

object

Structural results obtained using the Craig-Bampton order reduction method,
specified as a `ReducedStructuralModel`

object.

`u`

— Displacement

matrix

Displacement, specified as a matrix. The number of rows in the matrix must equal the
sum of the numbers of interface degrees of freedom and the number of modes. The
*x*-displacements at the retained degrees of freedom must appear
first, then the *y*-displacements, and, for a 3-D geometry,
*z*-displacements, followed by the generalized modal degrees of
freedom. The number of columns must equal the number of elements in
`tlist`

.

**Data Types: **`double`

`ut`

— Velocity

matrix

Velocity, specified as a matrix. The number of rows in the matrix must equal the sum
of the numbers of interface degrees of freedom and the number of modes. The
*x*-velocities at the retained degrees of freedom must appear first,
then the *y*-velocities, and, for a 3-D geometry,
*z*-velocities, followed by the generalized modal degrees of freedom.
The number of columns must equal the number of elements in
`tlist`

.

**Data Types: **`double`

`utt`

— Acceleration

matrix

Acceleration, specified as a matrix. The number of rows in the matrix must equal the
sum of the numbers of interface degrees of freedom and the number of modes. The
*x*-accelerations at the retained degrees of freedom must appear
first, then the *y*-accelerations, and, for a 3-D geometry,
*z*-accelerations, followed by the generalized modal degrees of
freedom. The number of columns must equal the number of elements in
`tlist`

.

**Data Types: **`double`

`tlist`

— Solution times for solving reduced-order model

real vector

Solution times for solving the reduced-order model, specified as a real vector.

**Data Types: **`double`

`Rtherm`

— Reduced-order thermal model

`ReducedThermalModel`

object

Reduced-order thermal model, specified as a `ReducedThermalModel`

object.

`u_therm`

— Temperature in modal coordinates

matrix

Temperature in modal coordinates, specified as a matrix. The number of rows in the
matrix must equal the number of modes. The number of columns must equal the number of
elements in `tlist`

.

**Data Types: **`double`

## Output Arguments

`structuralresults`

— Transient structural analysis results

`TransientStructuralResults`

object

Transient structural analysis results, returned as a `TransientStructuralResults`

object. The object contains the displacement,
velocity, and acceleration values at the nodes of the triangular or tetrahedral mesh
generated by `generateMesh`

.

`thermalresults`

— Transient thermal analysis results

`TransientThermalResults`

object

Transient thermal analysis results, returned as a `TransientThermalResults`

object. The object contains the temperature and
gradient values at the nodes of the triangular or tetrahedral mesh generated by
`generateMesh`

.

## Version History

**Introduced in R2019b**

### R2022a: ROM support for thermal analysis

`reconstructSolution`

now also reconstructs transient thermal
solutions.

## MATLAB Command

You clicked a link that corresponds to this MATLAB command:

Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)