Welcome to My Mechanical Engineering Site

This site has articles about stress analysis and other aspects of mechanical CAE, as well as my software stressRefine. StressRefine automatically extracts more accurate stresses from finite element models without remeshing (see the software page).

I am a mechanical engineer, Ph.D. Stanford, 1980, and was co-inventor and one of the lead developers of Mechanica, a pioneering easy to use p-adaptive finite element code. More details on my background are on the about page.

Coming Soon: Open-Source for stressRefine and Free P-Fea Course

This image has an empty alt attribute; its file name is image-19.png
Take the new free P-Fea course to learn what’s “under the hood”

Open-Source stressRefine

I am busy making the source code for stressRefine available opensource. It is free for academic, research, and engineering use per the GNU open-source General Public License (https://www.gnu.org/licenses/licenses.en.html#GPL). The source will be available in several versions:

  • Simple Version: does analysis of full FEA meshes only. This version does not do “breakout” (local region analysis). The elements and solver bookkeeping are simplified to be more easily readable, and are not implemented in parallel. This version is less efficient, but is recommended for starting out to understand the concepts
  • Intermediate Versions: Same as simple, but the elements and solver bookkeeping are optimized and implemented in parallel (multi-threaded using openMP).
  • Full Version: Same as intermediate version but also does “breakout” (local region analysis).

The above 3 versions are in the form of a single project that builds an executable. They are presented in Microsoft Visual Studio C++ projects, using the intel MKL library and the sparse solver pardiso included with it. I am also working on making versions available with Linux.

These versions are for users who want to use and customize a pre-existing code, or get stressRefine to run on a different platform. But I’d also like to support users who want to work p-adaptivity into your own codes. To this end, I’ll also be making the p-adaptive support routines available in a library. Be calling this library for support, you can easily convert existing conventional isoparametric elements to p elements. It is my hope that the library will alleviate the need to “reinvent the wheel” to create p-elements. A lot of the grunt work, for example of calculating 3D polynomial basis functions and assuring continuity across element and face boundaries is automatically done by the library. Students, developers and researchers can concentrate on more fun aspects like new error estimation techniques or novel applications.

Free Course

I’m also introducing a new free course. If you look at the menu at the top you’ll see it now has “course: P-Fea” . The course is intended to teach students and developers who want to learn how to write a p-adaptive code using the stressRefine library. It will start with the theory behind p-adaptivity, then take you through the practical aspects of creating p-adaptive elements using the library, and writing analysis routines that do p-adaptive analysis. It will continue on to show how to do analysis of local regions. This will all be presented first with simplified elements and solution assembly routines, then I’ll cover some code optimization and parallelization.

All three full executable versions will be available for upload by Friday, along with the first course lesson. The simple version in library form will be available by the end of next week.

Enhanced StressRefine Is Now Available

An new version of stressRefine is now available free for download. This has several enhancements over the previous version.

Nastran results can now be accessed in three different forms. Only Ascii (.f06) results are available in the previous version. Nastran binary (.op2) results can be accessed in the new version. In addition, data tables output by a postprocessor for the local region of interest can also be read, as described in my previous post. Using binary results is convenient because it is the default Nastran output, while outputting a data table is useful for large models because it allows data to be specified only for the region of interest, making translation faster.

StressRefine also now outputs stresses in Nastran ascii format, which can be displayed in your own postprocessor, in addition to the fringe plots available in the opensource postprocessor cgx which is bundled with stressRefine.

My favorite enhancement is automatic “breakout” of parts from an assembly for use in local region analysis which was described here. This is straightforward if you are working with a CAD program and have access to the model tree, but was tricky to implement when only the model data in a Nastran input file is available.

StressRefine Innovations

StressRefine has introduced multiple innovations which include:

  • Robust Error Estimation: StressRefine uses a blend of the Zienkiewicz-Zhu error norm [1] and traction discontinuities across element faces for error estimation. This is accurate for a wide range of polynomial order in elements. For models with localized stress concentrations, the polynomial order stays low in most of the model and is only raised locally where needed.
  • polynomial basis functions that are compatible with conventional isoparametric finite elements (such as 10 noded tets). StressRefine p elements can be mixed with conventional elements in the same mesh, and results and boundary can readily be exchanged between conventional meshes and stressRefine. These basis functions are packaged into a convenient library which allows any developer to use p-adaptive elements without “reinventing the wheel”.
  • Automatic detection and isolation of a variety of common singularities:
  • Automatic smoothing of nodal loads on surfaces: for example, if you apply a pressure load in a preprocessor, it may not show up as a pressure load in the Nastran input file, but instead as an “equivalent nodal load”. But it is only equivalent for quadratic elements. StressRefine has to notice that the all the nodes of a surface are loaded so this is really a distributed load, not multiple concentrated loads, to get a smooth result. An “energy-smoothing” procedure is then applied to recover the user’s original intent.
  • Automatic breakout modeling, including automatic extraction of parts from large assemblies.

StressRefine’s architecture is designed to easily allow it’s polynomial basis functions and elements to be incorporated into other FEA codes. After I do a bit of cleanup and documentation, the source code for stressRefine will shortly be made available.


  1. Zienkiewicz, O, and Zhu, J, “Superconvergent Patch Recovery and a posteriori error estimation in the finite element method”, Int. J. Num. Meth. Engr, 33, 1331, 1992

Using Data Tables From an FEA Postprocessor for Breakout Modeling

A data table is a easy way for a user to specify a region of a model from which to output desired quantities from a postprocessor to a file. I’m going to demonstrate how this is done for FEMAP, doing it with other postprocessors should be a similar process. For the purposes of breakout modeling the quantities needed to be output are node ids and displacements.

The motivation for doing this is that for large models, from which you want to analyze a relatively small local region, even the translation step can become time consuming. You may have a million element model and want more accurate results in a region of 10000 elements. It is inefficient to have an analysis tool like stressRefine wade through the results of a million elements to extract those needed in the local region.

A data table allows you to simply “rubber-band” select a region of interest. I’m going to show doing this with box-selecting.

So here’s how it goes in Femap:

Now we need to tell the Femap selector to box select nodes:

Note there are lots of other choices besides box. Now select a box a bit larger than the local region of interest. It doesn’t have to be precise, but it helps to include the whole part for “breakout by part”:

The data for all the nodes in the box will go into the data table. Now we specify what quantities to output (node ids and displacements).

Table is now ready to be saved:

This will all be in the stressRefine user’s guide that will be updated shortly. This procedure is a lot simpler to do than it is to explain, it takes just a few seconds once you have tried it a couple of times.


The example of breakout by part in an assembly from my previous post used a data table. Extraction of the breakout model took 2.4 secs using the data table, while it took over 30 additional seconds to translate the results and extract the model directly from the Nastran binary output results. This is for a model with less than 300,000 elements. The time savings will be much more for larger models.

Automatically Determining The Point of Maximum Stress In The Region Selected

In stressRefine, the user can specify the center of the local region using a node id. But this is not necessary by default:

The point of maximum stress in the region enclosed by the data table is used by default. Determining this point quickly without having to translate the results for the full model requires stressRefine to recreate the stresses in the elements that own the nodes in the data table. This is straightforward because displacements are known, and it is known that Nastran used quadratic elements (e.g. 10 node tets) in the solution. Calculating the stresses in the local region in the example above took much less than a second (for 55000+ nodes in the data table, owned by 33000+ nodes).

Breakout Analysis of Parts in An Assembly

I’ve recently discussed how local regions can be broken out of large models to more quickly calculate accurate stresses. These regions can have ragged boundaries as shown in the example below:

This does not affect the accuracy of the results as long as elements on the ragged boundary are treated as sacrificial. But sometimes the region of interest for more accurate stresses is a part in an assembly. For this case it is nice to have the entire part be the breakout region and figure out how to automatically extract it. Working with a CAD model, with access to the model tree makes this straightforward to implement. StressRefine’s interface with Solidworks Simulation does that, for example, by determining which surfaces form the boundary of a part using information from the Application Programming Interface (API). Through the API it is also easy to determine which surfaces are free boundaries and which are interfaces with other parts, from which the breakout boundary conditions can be determined. Here is an example using the Solidworks interface:

Working with Nastran models it appears the equivalent information is not available. However, recently I figured out how to use interface constraints (“bsurfs”) to cleanly breakout out parts from assemblies for use as local regions. Working from the center of the region of interest (which is either specified by the user by selecting a node, or is the point of max stress in the model), the mesh is traversed topologically, but does not go past faces with bsurfs.

The traversal to identify the part is shown progressing in the figure below. I also made a slide show that is a time-lapse sequence of the traversal which is shown here (you need to download it and run it with powerpoint or the powerpoint reader for it to show automatically, or in google slides hit view/present (ctrl-F5) to show the show but you’ll have to use page down or right arrow to advance it, it powerpoint it is automatic).

The overall model has 297,000 elements but extraction of the local region only took 2.4 secs, and stressRefine achieved a accurate p-adaptive solution of the local region in 18 secs. It helps that the local region is only 9963 elements, that’s the beauty of breakout analysis.

I think this “breakout by part” procedure is a useful enhancement to this technology. Just specify a point in the region of interest (or specify using the point of max stress in the model), and the appropriate part in an assembly will be automatically extracted as a local region. If that part has a very large mesh, stressRefine can also use a subset of the mesh in the part as the local region, again extracted automatically.

StressRefine and StressCheck

There are currently two well-known p-adaptive commercial FEA codes on the market, Mechanica and StressCheck.

Mechanica was an ambitious undertaking that those of us at Rasna and those who joined us later working for PTC worked hard on for 25 years. Towards the end we had added some pretty advanced features like large strain plasticity and large displacement sliding contact. We took pride in trying to make everything we did adaptive, by default. So all the underlying solutions in nonlinear analysis or optimization runs were p-adaptive. In addition, step size control in time-dependent solutions and nonlinear analysis were always adaptive. This did not always work, but it’s what we aimed for. Unfortunately PTC seems to be going in other directions and is not emphasizing Mechanica much these days. Recently they made a corporate deal with Ansys and our now selling Discovery Live under the name “Simulation Live”.

Another major commercial code on the market is StressCheck from ESRD. ESRD was founded by Professor Barna Szabo, who together with Professor Ivo Babuska invented the P-method, proved it converges, and worked out the important details of its implementation. ESRD recently had its 30th anniversary, and StressCheck is an excellent product which its customers use to achieve high quality results from FEA analyses. It has it’s own user interface and mesher. The user imports the CAD geometry and creates a high quality model in StressCheck, which then uses p-refinement to get accurate results. Users can also do breakout modeling. A local region, part of a large global model solved in another FEA code, is recreated in StressCheck. StressCheck provides tools to apply the boundary conditions to the local region from the solution to the global FEA analysis. StressCheck can do linear and various types of nonlinear analyses. I really like their CAE handbook, and they have a superb solution for composite joints.

In contrast, StressRefine is a niche product created for the sole purpose of being an accessory to other FEA codes, to extract more accurate stress results. It does not require the recreation of the model, instead working directly on an existing FEA model and using p-adaptivity to extract more accurate stresses from it. For large models it also does breakout analysis, but automates the process of creating the breakout region, applying interface conditions to it from the global model, and solving it adaptively for more accurate stress. So in it’s niche it requires less user intervention.

I would heartily recommend to anyone to use StressCheck. And perhaps consider my little niche product as an additional useful add-on.

Small Features and Stress Analysis

Small features are often ignored in finite element analysis, either by defeaturing them or not meshing finely in their vicinity. This can be fine for overall results like displacement or vibration where they often have little effect. But this may not the case for stress analysis. Consider this model with a slotted hole. The curves at the end of the slot do not become less important for stress as the slot becomes thinner. The stress trends upward with smaller radius, and in fact eventually becomes equivalent to a crack, with infinite stress, as the radius goes towards zero.

A famous real-world example of this it that some early versions of the DeHavilland Comet, the first commercial jetliner and a beautiful aircraft, crashed due to fatigue cracks emanating from the corners of their windows. The cornenrs had fillets that were too small, and the stress riser from that “small feature” caused the fatigue crack.

I was reminded of this recently while trying to modify automeshes in Femap. I came across this setting on their mesh sizing dialog:

If you miss this setting you can play around all you want with “Curvature-based mesh refinement”, “Surface Interior Mesh Growth”, etc and wonder why it is having no effect on the mesh near an important small feature like a fillet. You have to measure the size of the feature. If it is less than 0.0106 in the example above, it will be ignored, so you have to override this setting with a number smaller than your feature. I’m not trying to pick on Femap with this example, I’m sure there is a similar setting in other pre-processors. As mentioned above, the code doesn’t know your are interested in accurate stresses, and the feature may not have a big effect on other results.

What I would like to see in finite element structural analysis interfaces is a master setting like “stress results are important from this analysis”. Various default could key off this, such as turning “midside nodes on surfaces” on and “Curvature-based mesh refinement” . And the default setting of “max size of small feature” could be made smaller.

And I think this is an issue that users who need good stress results should be aware of. I’ve meshed large models and the mesh appears to be of good quality until you zoom in near important features. It may not be following the geometry well there, and adjacent elements may be large compared to the size of the feature.

In the example below, the overall mesh looks good until we zoom in near the fillet around the boss, and the mesh is way too coarse there. In the lower right zoomed in image, telling the mesher that “max size of small feature” should be less than 1 mm (the radius of the fillet) and turning on Curvature-based mesh refinement were needed to get this better mesh.