Experiences with Developing a “Somewhat Large” ACT Extension in ANSYS

With each release of ANSYS the customization toolkit continues to evolve and grow.  Recently I developed what I would categorize as a decent sized ACT extension.    My purpose in this post is to highlight a few of the techniques and best practices that I learned along the way.

Why I chose C#?

Most ACT extensions are written in Python.  Python is a wonderfully useful language for quickly prototyping and building applications, frankly of all shapes and sizes.  Its weaker type system, plethora of libraries, large ecosystem and native support directly within the ACT console make it a natural choice for most ACT work.  So, why choose to move to C#?

The primary reasons I chose to use C# instead of python for my ACT work were the following:

  1. I prefer the slightly stronger type safety afforded by the more strongly typed language. Having a definitive compilation step forces me to show my code first to a compiler.  Only if and when the compiler can generate an assembly for my source do I get to move to the next step of trying to run/debug.  Bugs caught at compile time are the cheapest and generally easiest bugs to fix.  And, by definition, they are the most likely to be fixed.  (You’re stuck until you do…)
  2. The C# development experience is deeply integrated into the Visual Studio developer tool. This affords not only a great editor in which to write the code, but more importantly perhaps the world’s best debugger to figure out when and how things went wrong.   While it is possible to both edit and debug python code in Visual Studio, the C# experience is vastly superior.

The Cost of Doing ACT Business in C#

Unfortunately, writing an ACT extension in C# does incur some development cost in terms setting up the development environment to support the work.  When writing an extension solely in Python you really only need a decent text editor.  Once you setup your ACT extension according to the documented directory structure protocol, you can just edit the python script files directly within that directory structure.  If you recall, ACT requires an XML file to define the extension and then a directory with the same name that contains all of the assets defining the extension like scripts, images, etc…  This “defines” the extension.

When it comes to laying out the requisite ACT extension directory structure on disk, C# complicates things a bit.  As mentioned earlier, C# involves a compilation step that produces a DLL.  This DLL must then somehow be loaded into Mechanical to be used within the extension.  To complicate things a little further, Visual Studio uses a predefined project directory structure that places the build products (DLLs, etc…) within specific directories of the project depending on what type of build you are performing.   Therefore the compiled DLL may end up in any number of different directories depending on how you decide to build the project.  Finally, I have found that the debugging experience within Visual Studio is best served by leaving the DLL located precisely wherever Visual Studio created it.

Here is a summary list of the requirements/problems I encountered when building an ACT extension using C#

  1. I need to somehow load the produced DLL into Mechanical so my extension can use it.
  2. The DLL that is produced during compilation may end up in any number of different directories on disk.
  3. An ACT Extension must conform to a predefined structural layout on the filesystem. This layout does not map cleanly to the Visual studio project layout.
  4. The debugging experience in Visual Studio is best served by leaving the produced DLL exactly where Visual Studio left it.

The solution that I came up with to solve these problems was twofold.

First, the issue of loading the proper DLL into Mechanical was solved by using a combination of environment variables on my development machine in conjunction with some Python programming within the ACT main python script.  Yes, even though the bulk of the extension is written in C#, there is still a python script to sort of boot-load the extension into Mechanical.  More on that below.

Second, I decided to completely rebuild the ACT extension directory structure on my local filesystem every time I built the project in C#.  To accomplish this, I created in visual studio what are known as post-build events that allow you to specify an action to occur automatically after the project is successfully built.  This action can be quite generic.  In my case, the “action” was to locally run a python script and provide it with a few arguments on the command line.  More on that below.

Loading the Proper DLL into Mechanical

As I mentioned above, even an ACT extension written in C# requires a bit of Python code to bootstrap it into Mechanical.  It is within this bit of Python that I chose to tackle the problem of deciding which dll to actually load.  The code I came up with looks like the following:

Essentially what I am doing above is querying for the presence of a particular environment variable that is on my machine.  (The assumption is that it wouldn’t randomly show up on end user’s machine…) If that variable is found and its value is 1, then I determine whether or not to load a debug or release version of the DLL depending on the type of build.  I use two additional environment variables to specify where the debug and release directories for my Visual Studio project exist.  Finally, if I determine that I’m running on a user’s machine, I simply look for the DLL in the proper location within the extension directory.  Setting up my python script in this way enables me to forget about having to edit it once I’m ready to share my extension with someone else.  It just works.

Rebuilding the ACT Extension Directory Structure

The final piece of the puzzle involves rebuilding the ACT extension directory structure upon the completion of a successful build.  I do this for a few different reasons.

  1. I always want to have a pristine copy of my extension laid out on disk in a manner that could be easily shared with others.
  2. I like to store all of the various extension assets, like images, XML files, python files, etc… within the Visual Studio Project. In this way, I can force the project to be out of date and in need of a rebuild if any of these files change.  I find this particularly useful for working with the XML definition file for the extension.
  3. Having all of these files within the Visual Studio Project makes tracking thing within a version control system like SVN or git much easier.

As I mentioned before, to accomplish this task I use a combination of local python scripting and post build events in Visual Studio.  I won’t show the entire python code, but essentially what it does is programmatically work through my local file system where the C# code is built and extract all of the files needed to form the ACT extension.  It then deletes any old extension files that might exist from a previous build and lays down a completely new ACT extension directory structure in the specified location.  The definition of the post build event is specified within the project settings in Visual Studio as follows:

As you can see, all I do is call out to the system python interpreter and pass it a script with some arguments.  Visual Studio provides a great number of predefined variables that you can use to build up the command line for your script.  So, for example, I pass in a string that specifies what type of build I am currently performing, either “Debug” or “Release”.  Other strings are passed in to represent directories, etc…

The Synergies of Using Both Approaches

Finally, I will conclude with a note on the synergies you can achieve by using both of the approaches mentioned above.  One of the final enhancements I made to my post build script was to allow it to “edit” some of the text based assets that are used to define the ACT extension.  A text based asset is something like an XML file or python script.  What I came to realize is that certain aspects of the XML file that define the extension need to be different depending upon whether or not I wish to debug the extension locally or release the extension for an end user to consume.  Since I didn’t want to have to remember to make those modifications before I “released” the extension for someone else to use, I decided to encode those modifications into my post build script.  If the post build script was run after a “debug” build, I coded it to configure the extension for optimal debugging on my local machine.  However, if I built a “release” version of the extension, the post build script would slightly alter the XML definition file and the main python file to make it more suitable for running on an end user machine.   By automating it in this way, I could easily build for either scenario and confidently know that the resulting extension would be optimally configured for the particular end use.

Conclusions

Now that I have some experience in writing ACT extensions in C# I must honestly say that I prefer it over Python.  Much of the “extra plumbing” that one must invest in in order to get a C# extension up and running can be automated using the techniques described within this post.  After the requisite automation is setup, the development process is really straightforward.  From that point onward, the increased debugging fidelity, added type safety and familiarity a C based language make the development experience that much better!  Also, there are some cool things you can do in C# that I’m not 100% sure you can accomplish in Python alone.  More on that in later posts!

If you have ideas for an ACT extension to better serve your business needs and would like to speak with someone who has developed some extensions, please drop us a line.  We’d be happy to help out however we can!

 

Reading ANSYS Mechanical Result Files (.RST) from C/C++ (Part 3)

ansys-fortran-to-c-cpp-1-00In the last post of this series I illustrated how I handled the nested call structure of the procedural interface to ANSYS’ BinLib routines.  If you recall, any time you need to extract some information from an ANSYS result file, you have to bracket the function call that extracts the information with a *Begin and *End set of function calls.  These two functions setup and tear down internal structures needed by the FORTRAN library.  I showed how I used RAII principles in C++ along with a stack data structure to manage this pairing implicitly.  However, I have yet to actually read anything truly useful off of the result file!  This post centers on the design of a set of C++ iterators that are responsible for actually reading data off of the file.  By taking the time to write iterators, we expose the ANSYS RST library to many of the algorithms available within the standard template library (STL), and we also make our own job of writing custom algorithms that consume result file data much easier.  So, I think the investment is worthwhile.

If you’ve programmed in C++ within the last 10 years, you’ve undoubtedly been exposed to the standard template library.  The design of the library is really rather profound.  This image represents the high level design of the library in a pictorial fashion:

ansys-fortran-to-c-cpp-3-01

On one hand, the library provides a set of generic container objects that provide a robust implementation of many of the classic data structures available within the field of computer science.  The collection of containers includes things like arbitrarily sized contiguous arrays (vectors), linked lists, associative arrays, which are implemented as either binary trees or as a hash container, as well as many more.  The set of containers alone make the STL quite valuable for most programmers.

On the other hand, the library provides a set of generic algorithms that encompass a whole suite of functionality defined in classic computer science.  Sorting, searching, rearranging, merging, etc… are just a handful of the algorithms provided by the library.  Furthermore, extreme care has been taken within the implementation of these algorithms such that an average programmer would hard pressed to produce something safer and more efficient on their own.

However, the real gem of the standard library are iterators.  Iterators bridge the gap between generic containers on one side and the generic algorithms on the other side.  Need to sort a vector of integers, or a double ended queue of strings?  If so, you just call the same sort function and pass it a set of iterators.  These iterators “know” how to traverse their parent container.  (Remember containers are the data structures.)

So, what if we could write a series of iterators to access data from within an ANSYS result file?  What would that buy us?  Well, depending upon which concepts our iterators model, having them available would open up access to at least some of the STL suite of algorithms.  That’s good.  Furthermore, having iterators defined would open up the possibility of providing range objects.  If we can provide range objects, then all of the sudden range based for loops are possible.  These types of loops are more than just syntactic sugar.  By encapsulating the bounds of iteration within a range, and by using iterators in general to perform the iteration, the burden of a correct implementation is placed on the iterators themselves.  If you spend the time to get the iterator implementation correct, then any loop you write after that using either the iterators or better yet the range object will implicitly be correct from a loop construct standpoint.  Range based for loops also make your code cleaner and easier to reason about locally.

Now for the downside…  Iterators are kind of hard to write.  The price for the flexibility they offer is paid for in the amount of code it takes to write them.  Again, though, the idea is that you (or, better yet somebody else) writes these iterators once and then you have them available to use from that point forward.

Because of their flexibility, standard conformant iterators come in a number of different flavors.  In fact, they are very much like an ice cream Sunday where you can pick and choose what features to mix in or add on top.  This is great, but again it makes implementing them a bit of a chore.  Here are some of the design decisions you have to answer when implementing an iterator:

Decision Options Choice for RST Reader Iterators
Dereference Data Type Anything you want Special structs for each type of iterator.
Iteration Category 1.       Forward iterator
2.       Single pass iterator
3.       Bidirectional iterator
4.       Random access iterator
Forward, Single Pass

Iterators syntactically function much like pointers in C or C++.  That is, like a pointer you can increment an iterator with the ++ operator, you can dereference an iterator with the * operator and you can compare two iterators for equality.  We will talk more about incrementing and comparisons in a minute, but first let’s focus on dereferencing.  One thing we have to decide is what type of data the client of our iterator will receive when they dereference it.  My choice is to return a simple C++ structure with data members for each piece of data.  For example, when we are iterating over the nodal geometry, the RST file contains the node number, the nodal coordinates and the nodal rotations.  To represent this data, I create a structure like this:ansys-fortran-to-c-cpp-3-02

I think this is pretty self-explanatory.  Likewise, if we are iterating over the element geometry section of an RST file, there is quite a bit of useful information for each element.  The structure I use in that case looks like this:

ansys-fortran-to-c-cpp-3-03

 

Again, pretty self-explanatory.  So, when I’m building a node geometry iterator, I’m going to choose the NodalCoordinateData structure as my dereference type.

The next decision I have to make is what “kind” of iterator I’m going to create.  That is, what types of “iteration” will it support?  The C++ standard supports a variety of iterator categories.  You may be wondering why anyone would ever care about an “iteration category”?  Well, the reason is fundamental to the design of the STL.   Remember that the primary reason iterators exist is to provide a bridge between generic containers and generic algorithms.  However, any one particular algorithm may impose certain requirements on the underlying iterator for the algorithm to function appropriately.

Take the algorithm “sort” for example.  There are, in fact, lots of different “sort” algorithms.  The most efficient versions of the “sort” algorithm require that an iterator be able to jump around randomly in constant time within the container.  If the iterator supports jumping around (a.k.a. random access) then you can use it within the most efficient sort algorithm.   However, there are certain kinds of iterators that don’t support jumping around.  Take a linked list container as an example.  You cannot randomly jump around in a linked list in constant time.  To get to item B from item A you have to follow the links, which means you have to jump from link to link to link, where each jump takes some amount of processing time.  This means, for example, there is no easy way to cut a linked list exactly in half even if you know how many items in total are in the list.  To cut it in half you have to start at the beginning and follow the links until you’ve followed size/2 number of links.  At that point you are at the “center” of the list.  However, with an array, you simply choose an index equal to size/2 and you automatically get to the center of the array in one step.  Many sort algorithms, as an example, obtain their efficiency by effectively chopping the container into two equal pieces and recursively sorting the two pieces.  You lose all that efficiency if you have to walk out to the center.

If you look at the “types” of iterators in the table above you will see that they build upon one another.  That is, at the lowest level, I have to answer the question, can I just move forward one step?  If I can’t even do that, then I’m not an iterator at all.  After that, assuming I can move forward one step, can I only go through the range once, or can I go over the range multiple times?  If I can only go over the range once, I’m a single pass iterator.  Truthfully, the forward iterator concept and the single pass iterator concept form level 1A and 1B of the iterator hierarchy.  The next higher level of functionality is a bidirectional iterator.  This type of iterator can go forward and backwards one step in constant time.  Think of a doubly linked list.  With forward and backward links, I can go either direction one step in constant time.  Finally, the most flexible iterator is the random access iterator.  These are iterators that really are like raw pointers.  With a pointer you can perform pointer arithmetic such that you can add an arbitrary offset to a base pointer and end up at some random location in a range.  It’s up to you to make sure that you stay within bounds.  Certain classes of iterators provide this level of functionality, namely those associated with vectors and deques.

So, the question is what type of iterator should we support?  Perusing through the FORTRAN code shipped with ANSYS, there doesn’t appear to be an inherent limitation within the functions themselves that would preclude random access.  But, my assumption is that the routines were designed to access the data sequentially.  (At least, if I were the author of the functions that is what I would expect clients to do.)  So, I don’t know how well they would be tested regarding jumping around.  Furthermore, disk controllers and disk subsystems are most likely going to buffer the data as it is read, and they very likely perform best if the data access is sequential.  So, even if it is possible to randomly jump around on the result file, I’m not sold on it being a good idea from a performance standpoint.  Lastly, I explicitly want to keep all of the data on the disk, and not buffer large chunks of it into RAM within my library.  So, I settled on expressing my iterators as single pass, forward iterators.  These are fairly restrictive in nature, but I think they will serve the purpose of reading data off of the file quite nicely.

Regarding my choice to not buffer the data, let me pause for a minute and explain why I want to keep the data on the disk. First, in order to buffer the data from disk into RAM you have to read the data off of the disk one time to fill the buffer.  So, that process automatically incurs one disk read.  Therefore, if you only ever need to iterate over the data once, perhaps to load it into a more specialized data structure, buffering it first into an intermediate RAM storage will actually slow you down, and consume more RAM.  The reason for this is that you would first iterate over the data stored on the disk and read it into an intermediate buffer.  Then, you would let your client know the data is ready and they would iterate over your internal buffer to access the data.  They might as well just read the data off the disk themselves! If the end user really wants to buffer the data, that’s totally fine.  They can choose to do that themselves, but they shouldn’t have to pay for it if they don’t need it.

Finally, we are ready to implement the iterators themselves.  To do this I rely on a very high quality open source library called Boost.  Boost has within it a library called iterator_facade that takes care of supplying most all of the boilerplate code needed to create a conformant iterator.  Using it has proven to be a real time saver.  To define the actual iterator, you derive your iterator class from iterator_facade and pass it a few template parameters.  One is the category defining the type of iterator you are creating.  Here is the definition for the nodal geometry iterator:

ansys-fortran-to-c-cpp-3-04

You can see that there are a few private functions here that actually do all of the work.  The function “increment” is responsible for moving the iterator forward one spot.  The function “equal” is responsible for determining if two different iterators are in fact equal.  And the function “dereference” is used to return the data associated with the current iteration spot.  You will also notice that I locally buffer the single piece of data associated with the current location in the iteration space inside the iterator.  This is stored in the pData member function.  I also locally store the current index.   Here are the implementations of the functions just mentioned:

ansys-fortran-to-c-cpp-3-05

First you can see that testing iterator equality is easy.  All we do is just look to see if the two iterators are pointing to the same index.  If they are, we define them as equal. (Note, an important nuance of this is that we don’t test to see if their buffered data is equal, just their index.  This is important later on.)  Likewise, increment is easy to understand as well.  We just increase the index by one, and then buffer the new data off of disk into our local storage.  Finally, dereference is easy too.  All we do here is just return a reference to the local data store that holds this one node’s data.  The only real work occurs in the readData() function.  Inside that function you will see the actual call to the CResRdNode() function.  We pass that function our current index and it populates an array of 6 doubles with data and returns the actual node number.  After we have that, we simply parse out of that array of 6 doubles the coordinates and rotations, storing them in our local storage.  That’s all there is to it.  A little more work, but not bad.

With these handful of operations, the boost iterator_facade class can actually build up a fully conformant iterator will all the proper operator overloads, etc… It just works.  Now that we have iterators, we need to provide a “begin” and “end” function just like the standard containers do.  These functions should return iterators that point to the beginning of our data set and to one past the end of our data set.  You may be thinking to yourself, wait a minute, how to we provide an “end” iterator without reading the whole set of nodes?  The reality is, we just need to provide an iterator that “equality tests” to be equal to the end of our iteration space?  What does that mean?  Well, what it means is that we just need to provide an iterator that, when compared to another iterator which has walked all the way to the end, it passes the “equal” test.  Look at the “equal” function above.  What do two iterators need to have in common to be considered equal?  They need to have the same index.  So, the “end” iterator simply has an index equal to one past the end of the iteration space.  We know how big our iteration space is because that is one of the pieces of metadata supplied by those ResRd*Begin functions.  So, here are our begin/end functions to give us a pair of conformant iterators.

ansys-fortran-to-c-cpp-3-06

Notice, that the nodes_end() function creates a NodeIterator and passes it an index that is one past the maximum number of nodes that have coordinate data stored on file.  You will also notice, that it does not have a second Boolean argument associated with it.  I use that second argument to determine if I should immediately read data off of the disk when the iterator is constructed.  For the begin iterator, I need to do that.  For the end iterator, I don’t actually need to cache any data.  In fact, no data for that node is defined on disk.  I just need a sentinel iterator that is one past the iteration space.

So, there you have it.  Iterators are defined that implicitly walk over the rst file pulling data off as needed and locally buffering one piece of it.  These iterators are standard conformant and thus can be used with any STL algorithm that accepts a single pass, read only, forward iterator.  They are efficient in time and storage.  There is, though, one last thing that would be nice.  That is to provide a range object so that we can have our cake and eat it too.  That is, so we can write these C++11 range based for loops.  Like this:ansys-fortran-to-c-cpp-3-07

To do that we need a bit of template magic.  Consider this template and type definition:

ansys-fortran-to-c-cpp-3-08

There is a bit of machinery that is going on here, but the concept is simple.  I just want the compiler to stamp out a new type that has a “begin()” and “end()” member function that actually call my “nodes_begin()” and “nodes_end()” functions.  That is what this template will do.  I can also create a type that will call my “elements_begin()” and “elements_end()” function.  Once I have those types, creating range objects suitable for C++11 range based for loops is a snap.  You just make a function like this:

ansys-fortran-to-c-cpp-3-09

 

This function creates one of these special range types and passes in a pointer to our RST reader.  When the compiler then sees this code:

ansys-fortran-to-c-cpp-3-10

It sees a range object as the return type of the “nodes()” function.  That range object is compatible with the semantics of range based for loops, and therefore the compiler happily generates code for this construction.

Now, after all of this work, the client of the RST reader library can open a result file, select something of interest, and start looping over the items in that category; all in three lines of code.  There is no need to understand the nuances of the binlib functions.  But best of all, there is no appreciable performance hit for this abstraction.  At the end of the day, we’re not computationally doing anything more than what a raw use of the binlib functions would perform.  But, we’re adding a great deal of type safety, and, in my opinion, readability to the code.  But, then again, I’m only slightly partial to C++.  Your mileage may vary…

Reading ANSYS Mechanical Result Files (.RST) from C/C++ (Part 2)

ansys-fortran-to-c-cpp-1-00In the last post in this series I illustrated how you can interface C code with FORTRAN code so that it is possible to use the ANSYS, Inc. BinLib routines to read data off of an ANSYS result file within a C or C++ program.  If you recall, the routines in BinLib are written in FORTRAN, and my solution was to use the interoperability features of the Intel FORTRAN compiler to create a shim library that sits between my C++ code and the BinLib code. In essence, I replicated all of the functions in the original BinLib library, but gave them a C interface. I call this library CBinLib.

You may remember from the last post that I wanted to add a more C++ friendly interface on top of the CBinLib functions.  In particular, I showed this piece of code, and alluded to an explanation of how I made this happen.  This post covers the first half of what it takes to make the code below possible.

ansys-fortran-to-c-cpp-2-01

What you see in the code above is the use of C++11 range based “for loops” to iterate over the nodes and elements stored on the result file.  To accomplish this we need to create conformant STL style iterators and ranges that iterate over some space.  I’ll describe the creation of those in a subsequent post.  In this post, however, we have to tackle a different problem.  The solution to the problem is hidden behind the “select” function call shown in the code above.  In order to provide some context for the problem, let me first show you the calling sequence for the procedural interface to BinLib.  This is how you would use BinLib if you were programming in FORTRAN or if you were directly using the CBinLib library described in the previous post.  Here is the recommended calling sequence:

ansys-fortran-to-c-cpp-2-02

You can see the design pattern clearly in this skeleton code.  You start by calling ResRdBegin, which gives you a bunch of useful data about the contents of the file in general.  Then, if you want to read geometric data, you need to call ResRdGeomBegin, which gives you a little more useful metadata.  After that, if you want to read the nodal coordinate data you call ResRdNodeBegin followed by a loop over nodes calling ResRdNode, which gives you the actual data about individual nodes, and then finally call ResRdNodeEnd.  If at that point you are done with reading geometric data, you then call ResRdGeomEnd.  And, if you are done with the file you call ResRdEnd.

Now, one thing jumps off the page immediately.  It looks like it is important to pair up the *Begin and*End calls.  In fact, if you peruse the ResRd.F FORTRAN file included with the ANSYS distribution you will see that in many of the *End functions, they release resources that were allocated and setup in the corresponding *Begin function.

So, if you forget to call the appropriate *End, you might leak resources.  And, if you forget to call the appropriate *Begin, things might not be setup properly for you to iterate.  Therefore, in either case, if you fail to call the right function, things are going to go badly for you.

This type of design pattern where you “construct” some scaffolding, do some work, and then “destruct” the scaffolding is ripe for abstracting away in a C++ type.  In fact, one of the design principles of C++ known as RAII (Resource Acquisition Is Initialization) maps directly to this problem.  Imagine that we create a class in which in the constructor of the class we call the corresponding *Begin function.  Likewise, in the destructor of the class we call the corresponding *End function.  Now, as long as we create an instance of the class before we begin iterating, and then hang onto that instance until we are done iterating, we will properly match up the *Begin, *End calls.  All we have to do is create classes for each of these function pairs and then create an instance of that class before we start iterating.  As long as that instance is still alive until we are finished iterating, we are good.

Ok, so abstracting the *Begin and *End function pairs away into classes is nice, but how does that actually help us?  You would still have to create an instance of the class, hold onto it while you are iterating, and then destroy it when you are done.  That sounds like more work than just calling the *Begin, *End directly.  Well, at first glance it is, but let’s see if we can use the paradigm more intelligently.  For the rest of this article, I’ll refer to these types of classes as BeginEnd classes, though I call them something different in the code.

First, what we really want is to fire and forget with these classes.  That is, we want to eliminate the need to manually manage the lifetime of these BeginEnd classes.  If I don’t accomplish this, then I’ve failed to reduce the complexity of the *Begin and *End requirements.  So, what I would like to do is to create the appropriate BeginEnd class when I’m ready to extract a certain type of data off of the file, and then later on have it magically delete itself (and thus call the appropriate *End function) at the right time.  Now, one more complication.  You will notice that these *Begin and*End function pairs are nested.  That is, I need to call ResRdGeomBegin before I call ResRdNodeBegin.  So, not only do I want a solution that allows me to fire and forget, but I want a solution that manages this nesting.

Whenever you see nesting, you should think of a stack data structure.  To increase the nesting level, you push an item onto the stack.  To decrease the nesting level, you pop and item off of the stack.  So, we’re going to maintain a stack of these BeginEnd classes.  As an added benefit, when we introduce a container into the design space, we’ve introduced something that will control object lifetime for us.  So, this stack is going to serve two functions for us.  It’s going to ensure we call the *Begin’s and *End’s in the right nested order, and second, it’s going to maintain the BeginEnd object lifetimes for us implicitly.

To show some code, here is the prototype for my pure virtual class that serves as a base class for all of the BeginEnd classes.  (In my code, I call these classes FileSection classes)

ansys-fortran-to-c-cpp-2-03

You can see that it is an interface class by noting the pure virtual function getLevel.  You will also notice that this function returns a ResultFileSectionLevel.  This is just an enum over file section types.  I like to use an enum as opposed to relying on RTTI.  Now, for each BeginEnd pair, I create an appropriate derived class from this base ResultFileSection class.  Within the destructor of each of the derived classes I call the appropriate *End function.  Finally, here is my stack data structure definition:

ansys-fortran-to-c-cpp-2-03p5

 

You can see that it is just a std::stack holding objects of the type SectionPtrT.  A SectionPtrT is a std::unique_ptr for objects convertible to my base section class.  This will enable the stack to hold polymorphic data, and the std::unique_ptr will manage the lifetime of the objects appropriately.   That is, when we pop and object off of the stack, the std::unique_ptr will make sure to call delete, which will call the destructor.  The destructor calls the appropriate *End function as we’ve mentioned before.

At this point, we’ve reduced our problem to managing items on a stack.  We’re getting closer to making our lives easier, trust me!  Let’s look at a couple of different functions to show how we pull these pieces together.  The first function is called reduceToLevelOrBegin(level).  See the code below:ansys-fortran-to-c-cpp-2-04

The operation of this function is fairly straightforward, yet it serves an integral role in our BeginEnd management infrastructure.   What this function does is it iteratively pops items off of our section stack until it either reaches the specified level, or it reaches the topmost ResRdBegin level.  Again, through the magic of C++ polymorphism, when an item gets popped off of the stack, eventually its destructor is called and that in turn calls the appropriate *End function.  So, what this function accomplishes is it puts us at a known level in the nested section hierarchy and, while doing so, ensures that any necessary *End functions are called appropriately to free resources on the FORTRAN side of things.  Notice that all of that happens automatically because of the type system in C++.  By popping items off of the stack, I implicitly clean up after myself.

The second function to consider is one of a family of similar functions.  We will look at the function that prepares the result file reader to read element geometry data off of the file.  Here it is:

ansys-fortran-to-c-cpp-2-05

You will notice that we start by reducing the nested level to either the “Geometry” level or the “Begin” level.  Effectively what this does is unwind any nesting we have done previously.  This is the machinery that makes “fire and forget” possible.  That is, whenever in ages past that we requested something to be read off of the result file, we would have pushed onto the stack a series of objects to represent the level needed to read the data in question.  Now that we wish to read something else, we unwind any previously existing nested Begin calls before doing so.   That is, we clean up after ourselves only when we ask to read a different set of data.  By waiting until we ask to read some new set of data to unwind the stack, we implicitly allow the lifetime of our BeginEnd classes to live beyond iteration.

At this point we have the stack in a known state; either it is at the Begin level or the Geometry level.  So, we simply call the appropriate *Begin functions depending on what level we are at, and push the appropriate type of BeginEnd objects onto the stack to record our traversal for later cleanup.  At this point, we are ready to begin iterating.  I’ll describe the process of creating iterators in the next post.  Clearly, there are lots of different select*** functions within my class.  I have chosen to make all of them private and have a single select function that takes an enum descriptor of what to select and some additional information for specifying result data.

One last thing to note with this design.  Closing the result file is easy. All that is required is that we simply fully unwind the stack.  That will ensure all of the appropriate FORTRAN cleanup code is called in the right order.  Here is the close function:

ansys-fortran-to-c-cpp-2-06

As you can see, cleanup is handled automatically.  So, to summarize, we use a stack of polymorphic data to manage the BeginEnd function calls that are necessary in the FORTRAN interface.  By doing this we ensure a level of safety in our class design.  Also, this moves us one step closer to this code:

ansys-fortran-to-c-cpp-2-07

In the next post I will show how I created iterators and range objects to enable clean for loops like the ones shown above.

Reading ANSYS Mechanical Result Files (.RST) from C/C++ (Part 1)

ansys-fortran-to-c-cpp-1-00Recently, I’ve encountered the need to read the contents of ANSYS Mechanical result files (e.g. file.rst, file.rth) into a C++ application that I am writing for a client. Obviously, these files are stored in a proprietary binary format owned by ANSYS, Inc.  Even if the format were published, it would be daunting to write a parser to handle it.  Fortunately, however, ANSYS supplies a series of routines that are stored in a library called BinLib which allow a programmer to access the contents of a result file in a procedural way.  That’s great!  But, the catch is the routines are written in FORTRAN… I’ve been programming for a long time now, and I’ll be honest, I can’t quite stomach FORTRAN.  Yes, the punch card days were before my time, but seriously, doesn’t a compiler have something better to do than gripe about what column I’m typing on… (Editor’s note: Matt does not understand the pure elegance of FORTRAN’s majestic simplicity. Any and all FORTRAN bashing is the personal opinion of Mr. Sutton and in no way reflects the opinion of PADT as a company or its owners. – EM) 

So, the problem shifts from how to read an ANSYS result file to how to interface between C/C++ and FORTRAN.  It turns out this is more complicated than it really should be, and that is almost exclusively because of the abomination known as CHARACTER(*) arrays.  Ah, FORTRAN… You see, if weren’t for the shoddy character of CHARACTER(*) arrays the mapping between the basic data types in FORTRAN and C would be virtually one for one. And thus, the mechanics of function calls could fairly easily be made to be identical between the two languages.   If the function call semantics were made identical, then sharing code between the two languages would be quite straightforward.  Alas, because a CHARACTER array has a kind of implicit length associated with it, the compiler has to do some kind of magic within any function signature that passes one or more of these arrays.  Some compilers hide parameters for the length and then tack them on to the end of the function call.  Some stuff the hidden parameters right after the CHARACTER array in the call sequence.  Some create a kind of structure that combines the length with the actual data into a special type. And then some compilers do who knows what…  The point is, there is no convention among FORTRAN compilers for handling the function call semantics, so there is no clean interoperability between C and FORTRAN.

Fortunately, the Intel FORTRAN compiler has created this markup language for FORTRAN that functions as an interoperability framework that enables FORTRAN to speak C and vice versa.  There are some limitations, however, which I won’t go into detail on here.  If you are interested you can read about it in the Intel FORTRAN compiler manual.  What I want to do is highlight an example of what this looks like and then describe how I used it to solve my problem.  First, an example:

ansys-fortran-to-c-cpp-1-01

What you see in this image is code for the first function you would call if you want to read an ANSYS result file.  There are a lot of arguments to this function, but in essence what you do is pass in the file name of the result file you wish to open (Fname), and if everything goes well, this function sends back a whole bunch of data about the contents of the file.  Now, this function represents code that I have written, but it is a mirror image of the ANSYS routine stored in the binlib library.

I’ve highlighted some aspects of the code that constitute part of the interoperability markup.  First of all, you’ll notice the markup BIND highlighted in red.  This markup for the FORTRAN function tells the compiler that I want it to generate code that can be called from C and I want the name of the C function to be “CResRdBegin”.  This is the first step towards making this function callable from C.  Next, you will see highlighted in blue something that looks like a comment.  This however, instructs the compiler to generate a stub in the exports library for this routine if you choose to compile it into a DLL.  You won’t get a .lib file when compiling this as a .dll without this attribute.  Finally, you see the ISO_C_BINDING and definition of the type of character data we can make interoperable.  That is, instead of a CHARACTER(261) data type, we use an array of single CHARACTER(1) data.  This more closely matches the layout of C, and allows the FORTRAN compiler to generate compatible code.  There is a catch here, though, and that is in the Title parameter.  ANSYS, Inc. defines this as an array of CHARACTER(80) data types.  Unfortunately, the interoperability stuff from Intel doesn’t support arrays of CHARACTER(*) data types.  So, we flatten it here into a single string.  More on that in a minute.

You will notice too, that there are markups like (c_int), etc… that the compiler uses to explicitly map the FORTRAN data type to a C data type.  This is just so that everything is explicitly called out, and there is no guesswork when it comes to calling the routine.  Now, consider this bit of code:

ansys-fortran-to-c-cpp-1-02

First, I direct your attention to the big red circle. Here you see that all I am doing is collecting up a bunch of arguments and passing them on to the ANSYS, Inc. routine stored in BinLib.lib.  You also should notice the naming convention.  My FORTRAN function is named CResRdBegin, whereas the ANSYS, Inc. function is named ResRdBegin.  I continue this pattern for all of the functions defined in the BinLib library.  So, this function is nothing more than a wrapper around the corresponding binlib routine, but it is annotated and constrained to be interoperable with the C programming language.  Once I compile this function with the FORTRAN compiler, the resulting code will be callable directly from C.

Now, there are a few more items that have to be straightened up.  I direct your attention to the black arrow.  Here, what I am doing is converting the passed in array of CHARACTER(1) data into a CHARACTER(*) data type.  This is because the ANSYS, Inc. version of this function expects that data type.  Also, the ANSYS, Inc. version needs to know the length of the file path string.  This is stored in the variable ncFname.  You can see that I compute this value using some intrinsics available within the compiler by searching for the C NULL character.  (Remember that all C strings are null terminated and the intent is to call this function from C and pass in a C string.)

Finally, after the call to the base function is made, the strings representing the JobName and Title must be converted back to a form compatible with C.  For the jobname, that is a fairly straightforward process.  The only thing to note is how I append the C_NULL_CHAR to the end of the string so that it is a properly terminated C string.

For the Title variable, I have to do something different.  Here I need to take the array of title strings and somehow represent that array as one string.  My choice is to delimit each title string with a newline character in the final output string.  So, there is a nested loop structure to build up the output string appropriately.

After all of this, we have a C function that we can call directly.  Here is a function prototype for this particular function.

ansys-fortran-to-c-cpp-1-03

So, with this technique in place, it’s just a matter of wrapping the remaining 50 functions in binlib appropriately!  Now, I was pleased with my return to the land of C, but I really wanted more.  The architecture of the binlib routines is quite easy to follow and very well laid out; however, it is very, very procedural for my tastes.  I’m writing my program in C++, so I would really like to hide as much of this procedural stuff as I can.   Let’s say I want to read the nodes and elements off of a result file.  Wouldn’t it be nice if my loops could look like this:

ansys-fortran-to-c-cpp-1-04

That is, I just do the following:

  1. Ask to open a file (First arrow)
  2. Tell the library I want to read nodal geometric data (Second arrow)
  3. Loop over all of the nodes on the RST file using C++11 range based for loops
  4. Repeat for elements

Isn’t that a lot cleaner?  What if we could do it without buffering data and have it compile down to something very close to the original FORTRAN code in size and speed?  Wouldn’t that be nice?  I’ll show you how I approached it in Part 2.

PID Thermostat Boundary Condition ACT Extension for ANSYS Mechanical

ANSYS-ACT-PID-ThermostatPADT is pleased to announce that we have uploaded a new ACT Extension to the ANSYS ACT App Store.  This new extension implements a PID based thermostat boundary condition that can be used within a transient thermal simulation.  This boundary condition is quite general purpose in nature.  For example, it can be setup to use any combination of (P)roportional (I)ntegral or (D)erivate control.   It supports locally monitoring the instantaneous temperature of any piece of geometry within the model.  For a piece of geometry that is associated with more than one node, such as an edge or a face, it uses a novel averaging scheme implemented using constraint equations so that the control law references a single temperature value regardless of the reference geometry.

ANSYS-ACT-PID-Thermostat-img1

The set-point value for the controller can be specified in one of two ways.  First, it can be specified as a simple table that is a function of time.  In this scenario, the PID ACT Extension will attempt to inject or remove energy from some location on the model such that a potentially different location of the model tracks the tabular values.   Alternatively, the PID thermostat boundary condition can be set up to “follow” the temperature value of a portion of the model.  This location again can be a vertex, edge or face and the ACT extension uses the same averaging scheme mentioned above for situations in which more than one node is associated with the reference geometry.  Finally, an offset value can be specified so that the set point temperature tracks a given location in the model with some nonzero offset.

ANSYS-ACT-PID-Thermostat-img2

For thermal models that require some notion of control the PID thermostat element can be used effectively.  Please do note, however, that the extension works best with the SI units system (m-kg-s).

Programming a Simple Polygon Editor

polygon-editor-icon-1Part of my job at PADT is writing custom software for our various clients.  We focus primarily on developing technical software for the engineering community, with a particular emphasis on tools that integrate with the ANSYS suite of simulation tools.  Frankly, writing software is my favorite thing to do at PADT, simply because software development is all about problem solving.

This morning I got to work on a fairly simple feature of a much larger tool that I am currently developing.  The feature I’m working on involves graphically editing polygons.  Why, you ask am I doing this?  Well, that I can’t say, but nonetheless I can share a particularly interesting problem (to me at least) that I got to take a swing at solving.  The problem is this:

When a user is editing a node in the polygon by dragging it around on the screen, how do you handle the case when they drop it on an existing node?

Consider this polygon I sketched out in a prototype of the tool.

polygon-editor-f01

What should happen if the user drags this node over on top of that node:polygon-editor-f02

Well, I think the most logical thing to do is that you merge the two nodes together.  Implementing that is pretty easy.  The slightly harder question is what to do with the remaining structure of the polygon?  For my use case, polygons have to be manifold in that no vertex is connected to more than two edges. (The polygons can be open and thus have two end vertices connected to only one edge.)  So, what part do you delete?  Well, my solution is that you delete the “smaller” part, where “smaller” is defined as the part that has the fewest nodes.  So, for example, this is what my polygon looks like after the “drop”

polygon-editor-f03Conceptually, this sounds pretty simple, but how do you do it programmatically?  To give some background, note that the nodes in my polygon class are stored in a simple, ordered C++ std::list<>.

Now, I use a std::list<> simply because I know I’m going to be inserting and deleting nodes at random places.  Linked lists are great for that, and for rendering, I have to walk the whole list anyway, so there’s no performance hit there.  Graphically, my data structure looks
something like this:polygon-editor-f04Pretty simple.  For a closed polygon, my class maintains a flag and simply draws an edge from the last node to the first node.

The rub comes when you start to realize that there are tons of different ways a user might try to merge nodes together in either an open or closed polygon.  I’ve illustrated a few below along with what nodes would need to be merged in the corresponding data structure.  In the data structure pictures, the red node is the target (the node on which the user will be dropping) and the green node is the one they are manipulating (the source node).

Here is one example:

polygon-editor-f05polygon-editor-f06

Here is another example:

polygon-editor-f07
polygon-editor-f09b

Finally, here is one more:polygon-editor-f08polygon-editor-f09

In all these examples, we have different “cases” that we need to handle.  For instance, in the first example the portion of the data structure we want to keep is the stuff between the source and target nodes.  So, the stuff on the “ends” of the list needs to be deleted.  In the middle case, we just need to merge the source and target together.  Finally, in the last case, the nodes between the source and target need to be deleted, whereas the stuff at the “ends” of the list need to be kept.

This simple type of problem causes shivers in many programmers, and I’ll admit, I was nervous at first that this problem was going to lead to a solution that handled each individual case respectively.  Nothing in all of programming is more hideous than that.  So, there has to be a simple way to figure out what part of the list to keep, and what part of the list to throw away.

Now, I’m sure this problem has been solved numerous times before, but I wanted to take a shot at it without googling.  (I still haven’t googled, yet… so if this is similar to any other approach, they get the credit and I just reinvented the wheel…)  I remember a long time ago listening to a C++ programmer espouse the wonders of the standard library’s algorithm section.  I vaguely remember him droning on about how wonderful the std::rotate algorithm is.  At the time, I didn’t see what all the fuss was about.  Now, I’m right there with him.  std::rotate is pretty awesome!

std::rotate is a simple algorithm.  Essentially what it does is it takes the first element in a list, pops it off the list and moves it to the rear of the list.  Everything else in the list shifts up one spot.  This is called a left rotate, because you can imagine the items in the list rotating to the left until they get to the front of the line, at which point they fall off and are put back on the end of the list.  (Using reverse iterators you can effectively perform a right rotate as well.)  So, how can we take advantage of this to simplify figuring out what needs to be deleted from our list of nodes?

Well, the answer is remarkably simple.  Once we locate the source and target nodes in the list, regardless of their relative position with respect to one another or to the ends of the list, we simply left rotate the list until the target becomes the head of the list.  That is, if we start with this:polygon-editor-f10We left rotate until we have this:polygon-editor-f11That’s great, but what does that buy us?  Well, now that one of the participating nodes is at the head of the list, our problem is much simpler because all of the nodes that we need to delete are now at either end of the list.  The only question left to answer is which end of the list do we trim off?  The answer to that question is trivial.  We simply trim off the shorter end of the list with respect to the source node (the green node in the diagram). The “lengths” of the two lists are defined as follows.  For the head section, it’s the number of nodes up to, but not including the source. (This section obviously includes the target node)  For the tail, it’s the number of nodes from the source to the end, including the source.  (This section includes the source node).  Since we define the two sections this way we are guaranteed to delete either the source or the target, but not both.  Its fine to delete either one of them, because at this point we’ve deemed the geometrically coincident, but we must not accidentally delete both!!

In the example just given, after the rotate, we would delete the head of the list.  However, let’s take a look at our first example.  Here is the original list:

polygon-editor-f12Here is the rotated list:polygon-editor-f13So, in this case, the “end” of the list (including the source) is the shortest.  If it is a tie, then it doesn’t matter, just pick one.  Interestingly enough, if the two nodes are adjacent in the original list, then the rotated list will look like either this:polygon-editor-f14 Or this, if the source is “before” the target in the original list:polygon-editor-f15In either case, the algorithm works unchanged, and we only delete one node.  It’s beautiful! (At least in my opinion…)  Modern C++ makes this type of code really clean and easy to write.  Here is the entire thing, including the search to located geometrically adjacent nodes as well as the merge. The standard library algorithms really help out!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
// Search lambda function for looking for any other node in the list that is
// coindicent to this node, except this node.
auto searchAdjacentFun = [this, pNode](const NodeListTool::AdjustNodePtrT &amp;pOtherNode)-&gt;bool
{
if (pNode-&gt;tag() == pOtherNode-&gt;tag()) return false;
return (QVector2D(pNode-&gt;pos() - pOtherNode-&gt;pos()).length() &lt; m_snapTolerance); }; auto targetLoc = std::find_if(m_nodes.begin(), m_nodes.end(), searchAdjacentFun); // If we don't find an adjacent node within the tolerance, then we can't merge if (targetLoc == m_nodes.end()) { return false; } // Tidy things up so that the source has exactly the same position as the target pNode-&gt;setPos((*targetLoc)-&gt;pos());
 
// Begin the merge by left rotating the target so that it is at the
// beginning of the list
std::rotate(m_nodes.begin(), targetLoc, m_nodes.end());
 
// Find this node in the list
auto searchThis = [this, pNode](const NodeListTool::AdjustNodePtrT &amp;pOtherNode)-&gt;bool
{
return (pNode-&gt;tag() == pOtherNode-&gt;tag());
};
auto sourceLoc = std::find_if(m_nodes.begin(), m_nodes.end(), searchThis);
 
// Now, figure out which nodes we are going to delete.
auto distToBeg = std::distance(m_nodes.begin(), sourceLoc);
auto distToEnd = std::distance(sourceLoc, m_nodes.end());
 
if (distToBeg &lt; distToEnd) { // If our source is closer to the beginning (which is the target) // than it is to the end of the list, then we need to delete // the nodes at the front of the list m_nodes.erase(m_nodes.begin(), sourceLoc); } else { // Otherwise, delete the nodes at the end of the list m_nodes.erase(sourceLoc, m_nodes.end()); } // Now, see if we still have more than 2 vertices if (m_nodes.size() &gt; 2) {
m_bClosed = true;
}
else {
m_bClosed = false;
}
return true;

Home Grown HPC on CUBE Systems

compute-cluster-1

A Little Project Background

Recently I’ve been working on developing a computer vision system for a long standing customer. We are developing software that enables them to use computers to “see” where a particular object is space, and accurately determine its precise location with respect to the camera. From that information, they can do all kinds of useful things.

In order to figure out where something is in 3D space from a 2D image you have to perform what is commonly referred to as pose estimation. It’s a highly interesting problem by itself, but it’s not something I want to focus on in detail here. If you are interested in obtaining more information, you can Google pose estimation or PnP problems. There are, however, a couple of aspects of that problem that do pertain to this blog article. First, pose estimation is typically a nonlinear, iterative process. (Not all algorithms are iterative, but the ones I’m using are.) Second, like any algorithm, its output is dependent upon its input; namely, the accuracy of its pose estimate is dependent upon the accuracy of the upstream image processing techniques. Whatever error happens upstream of this algorithm typically gets magnified as the algorithm processes the input.

The Problem I Wish to Solve

You might be wondering where we are going with HPC given all this talk about computer vision. It’s true that computer vision, especially image processing, is computationally intensive, but I’m not going to focus on that aspect. The problem I wanted to solve was this: Is there a particular kind of pattern that I can use as a target for the vision system such that the pose estimation is less sensitive to the input noise? In order to quantify “less sensitive” I needed to do some statistics. Statistics is almost-math, but just a hair shy. You can translate that statement as: My brain neither likes nor speaks statistics… (The probability of me not understanding statistical jargon is statistically significant. I took a p-test in a cup to figure that out…) At any rate, one thing that ALL statistics requires is a data set. A big data set. Making big data sets sounds like an HPC problem, and hence it was time to roll my own HPC.

The Toolbox and the Solution

My problem reduced down to a classic Monte Carlo type simulation. This particular type of problem maps very nicely onto a parallel processing paradigm known as Map-Reduce. The concept is shown below:
matt-hpc-1

The idea is pretty simple. You break the problem into chunks and you “Map” those chunks onto available processors. The processors do some work and then you “Reduce” the solution from each chunk into a single answer. This algorithm is recursive. That is, any single “Chunk” can itself become a new blue “Problem” that can be subdivided. As you can see, you can get explosive parallelism.

Now, there are tools that exist for this kind of thing. Hadoop is one such tool. I’m sure it is vastly superior to what I ended up using and implementing. However, I didn’t want to invest at this time in learning a specialized tool for this particular problem. I wanted to investigate a lower level tool on which this type of solution can be built. The tool I chose was node.js (www.nodejs.org).

I’m finding Node to be an awesome tool for hooking computers together in new and novel ways. It acts kind of like the post office in that you can send letters and messages and get letters and messages all while going about your normal day. It handles all of the coordinating and transporting. It basically sends out a helpful postman who taps you on the shoulder and says, “Hey, here’s a letter”. You are expected to do something (quickly) and maybe send back a letter to the original sender or someone else. More specifically, node turns everything that a computer can do into a “tap on the shoulder”, or an event. Things like: “Hey, go read this file for me.”, turns into, “OK. I’m happy to do that. I tell you what, I’ll tap you on the shoulder when I’m done. No need to wait for me.” So, now, instead of twiddling your thumbs while the computer spins up the harddrive, finds the file and reads it, you get to go do something else you need to do. As you can imagine, this is a really awesome way of doing things when stuff like network latency, hard drives spinning and little child processes that are doing useful work are all chewing up valuable time. Time that you could be using getting someone else started on some useful work. Also, like all children, these little helpful child processes that are doing real work never seem to take the same time to do the same task twice. However, simply being notified when they are done allows the coordinator to move on to other children. Think of a teacher in a class room. Everyone is doing work, but not at the same pace. Imagine if the teacher could only focus on one child at a time until that child fully finished. Nothing would ever get done!

Here is a little graph of our internal cluster at PADT cranking away on my Monte Carlo simulation.
matt-hpc-2

It’s probably impossible to read the axes, but that’s 1200+ cores cranking away. Now, here is the real kicker. All of the machines have an instance of node running on them, but one machine is coordinating the whole thing. The CPU on the master node barely nudges above idle. That is, this computer can manage and distribute all this work by barely lifting a finger.

Conclusion

There are a couple of things I want to draw your attention to as I wrap this up.

  1. CUBE systems aren’t only useful for CAE simulation HPC! They can be used for a wide range of HPC needs.
  2. PADT has a great deal of experience in software development both within the CAE ecosystem and outside of this ecosystem. This is one of the more enjoyable aspects of my job in particular.
  3. Learning new things is a blast and can have benefit in other aspects of life. Thinking about how to structure a problem as a series of events rather than a sequential series of steps has been very enlightening. In more ways than one, it is also why this blog article exists. My Monte Carlo simulator is running right now. I’m waiting on it to finish. My natural tendency is to busy wait. That is, spin brain cycles watching the CPU graph or the status counter tick down. However, in the time I’ve taken to write this article, my simulator has proceeded in parallel to my effort by eight steps. Each step represents generating and reducing a sample of 500,000,000 pose estimates! That is over 4 billion pose estimates in a little under an hour. I’ve managed to write 1,167 words…

CUBE_Logo_150w

ACT Extension for a PID Thermostat Controller (PART 2)

(View part one of this series here.)

So, I’ve done a little of this Workbench customization stuff in a past life. My customizations involved lots of Jscript, some APDL, sweat and tears. I literally would bang my head against Eric Miller’s office door jamb wondering (sometimes out loud) what I had done to be picked as the Workbench customization guy. Copious amounts of alcohol on weeknights helped some, but honestly it still royally sucked. Because of these early childhood scars, I’ve cringed at the thought of this ACT thingy until now. I figured I’d been there, done that and I had zero, and I mean zero desire to relive any of it.

So, I resisted ACT even after multiple “suggestions” from upper management that I figure out something to do with it. That was wrong of me; I should have been quicker to given ACT a fair shot. ACT is a whole bunch better than the bad ol’ days of JScript. How is it better? Well, it has documentation… Also, there are multiple helpful folks back at Canonsburg and elsewhere that know a few things about it. This is opposed to the days when just one (brilliant) guy in India named Rajiv Rath had ever done anything of consequence with JScript. (Without him, my previous customization efforts would simply have put me in the mad house. Oh, and he happens to know a thing or two about ACT as well…)

Look Ma! My First ACT Extension!

In this post we are going to rig up the PID thermostat boundary condition as a new boundary condition type in Mechanical. In ACT jargon, this is known as adding a pre-processing feature. I’m going to refer you primarily to the training material and documentation for ACT that you can obtain from the ANSYS website on their customer portal. I strongly suggest you log on to the portal and download the training material for version 15.0 now.

Planning the Extension

When we create an ACT extension we need to lay things out on the file system in a certain manner. At a high level, there are three categories of file system data that we will use to build our extension. These types of data are:

  1. Code. This will be comprised of Python scripts and in our case APDL scripts.
  2. XML. XML files are primarily used for plumbing and for making the rest of the world aware of who we are and what we do.
  3. Resources. These files are typically images, icons, etc… that spice things up a little bit.

Any single extension will use all three of these categories of files, and so organizing this mess becomes priority number one when building the extension. Fortunately, there is a convention that ACT imposes on us to organize the files. The following image depicts the structure graphically.

image

We will call our extension PIDThermostat. Therefore, everywhere ExtensionName appears in the image above, we will replace that with PIDThermostat in our file structure.

Creating a UI for our Load

The beauty of ACT is that it allows us to easily create professional looking user experiences for custom loads and results. Let’s start by creating a user interface for our load that looks like the following:

clip_image004

As you can see in the above user interface, all of the controls and inputs for our little PID controller that we designed in Part 1 of this blog series are present. Before we discuss how we create these user elements, let’s start with a description of what they each mean.

The first item in the UI is named Heat Source/Sink Location. This UI element presents to the user a choice between picking a piece of geometry from the model and specifying a named selection. Internal to the PID controller, this location represents where in the model we will attach the control elements such that they supply or remove energy from this location. ACT provides us two large areas of functionality here. First, it provides a way to graphically pick a geometric item from the model. Second, it provides routines to query the underlying mesh associated with this piece of geometry so that we can reference the nodes or elements in our APDL code. Each of these pieces of functionality is substantial in its size and complexity, yet we get it for free with ACT.

The second control is named Temperature Monitor Location. It functions similarly to the heat source/sink location. It gives the user the ability to specify where on the model the control element should monitor, or sense, the temperature. So, our PID controller will add or remove energy from the heat sink/source location to try to keep the monitor location at the specified set point temperature.

The third control group is named Thermostat Control Properties. This group aggregates the various parameters that control the functionality of the thermostat. That is, the user is allowed to specify gain values, and also what type of control is implemented.

The forth control group is named Thermostat Setpoint Properties. This group allows the user to specify how the setpoint changes with time. An interesting ACT feature is implemented with this control group. Based on the selection that the user makes in the “Setpoint Type” dropdown, a different control is presented below for the thermostat setpoint temperature. When the user selects, “User Specified Setpoint” then a control that provides a tabular input is presented. In this case, the user can directly input a table of temperature vs time data that specifies how the setpoint changes with time. However, if the user specifies “Use Model Entity as Setpoint” then the user is presented a geometry picker similar to the controls above to select a location in the model that will define the setpoint temperature. When this option is selected, the PID controller will function more like a follower element. That is, it will try to cause the monitor location to “follow” another location in the model by adding or removing energy from the heat source/sink location. The offset value allows the user to specify a DC offset that they would like to apply to the setpoint value. Internally, this offset term will be incorporated into the constraint equation averaging method to add in the DC offset.

Finally, the last control group allows the user to visualize plots of computed information regarding the PID controller after the solution is finished. Normally this would be presented in the results branch of the tree; however, the results I would like to present for these elements don’t map cleanly to the results objects in ACT. (At least, I can’t map them cleanly in my mind…) More detail on the results will be presented below.

So, now that we know what the control UI does, let’s look at how to specify it in ACT

Specifying the UI in XML

As mentioned at the beginning, ACT relies on XML to specify the UI for controls. The following XML snippet describes the entire UI.

<extension version=“1” name=“ThermalTools”>

  <guid shortid=“thermaltools”>852acb16-4731-4e91-bd00-b464be02b361</guid>

  <script src=“thermaltools.py” />

 

  <interface context=“Mechanical”>

    <images>images</images>

    <callbacks>

      <oninit>initThermalToolsToolbar</oninit>

    </callbacks>

    <toolbar name=“thermtools” caption=“Thermal Tools”>

      <entry name=“PIDThermostatLoad” icon=“ThermostatGray”

             caption=“PID Thermostat Control”>

        <callbacks>

          <onclick>createPIDThermostatLoad</onclick>

        </callbacks>

      </entry>

    </toolbar>

  </interface>

 

  <simdata context=“Mechanical”>

    <load name=“pidthermostat” version=“1” caption=“PID Thermostat”

    icon=“ThermostatWhite” issupport=“false” isload=“true” color=“#0000FF”>

      <callbacks>

        <getcommands location=“solve”>writePIDThermostatLoad</getcommands>

        <getcommands location=“post”>writePIDThermostatPost</getcommands>

        <onadd>addPIDThermostat</onadd>

        <onremove>removePIDThermostat</onremove>

        <onsuppress>suppressPIDThermostat</onsuppress>

        <onunsuppress>unsuppressPIDThermostat</onunsuppress>

        <canadd>canaddPIDThermostat</canadd>

        <oncleardata>cleardataPIDThermostat</oncleardata>

      </callbacks>

      <property name=“ConnectionGeo”  caption= “Heat Source/Sink Location”

                control=“scoping”>

        <attributes selection_filter=“vertex|edge|face” />

      </property>

      <property name=“MonitorGeo”  caption= “Temperature Monitor Location”

                control=“scoping”>

        <attributes selection_filter=“vertex|edge|face|body” />

      </property>

      <propertygroup name=“ControlProperties” caption=“Thermostat Control Properties”

                     display=“caption”>

        <property name=“ControlType” caption=“Control Type”

                  control=“select” default=“Both Heat Source and Sink”>

          <attributes options=“Heat Source,Heat Sink,Both Heat Source and Sink”/>

        </property>

        <property name=“ProportionalGain” caption=“Proportional Gain”

                  control=“float” default=“0”/>

        <property name=“IntegralGain” caption=“Integral Gain”

                  control=“float” default=“0”/>

        <property name=“DerivativeGain” caption=“Derivative Gain”

                  control=“float” default=“0”/>

      </propertygroup>

      <propertygroup name=“SetpointProperties” caption=“Thermostat Setpoint Properties”

                     display=“caption”>

        <propertygroup name=“SetpointType” display=“property” caption=“Setpoint Type”

                       control=“select” default=“User Specified Setpoint”>

          <attributes options=“User Specified Setpoint,Use Model Entity as Setpoint”/>

          <propertytable name=“SetPointTemp” caption=“Thermostat Set Point Temperature”

                         display=“worksheet” visibleon=“User Specified Setpoint”

                         control=“applycancel” class=“Worksheet.PropertyGroupEditor.PGEditor”>

            <property name=“Time” caption=“Time” unit=“Time” control=“float”></property>

            <property name=“SetPoint” caption=“Set Point Temperature”

                      unit=“Temperature” control=“float”></property>

          </propertytable>

          <property name=“SetpointGeo”  caption= “Setpoint Geometry”

                    visibleon=“Use Model Entity as Setpoint” control=“scoping”>

            <attributes selection_filter=“vertex|edge|face|body” />

          </property>

        </propertygroup>

        <property name=“SetpointOffset” caption=“Offset” control=“float” default=“0”/>

      </propertygroup>

      <propertygroup name=“Results” caption=“Thermostat Results” display=“caption”>

        <property name=“ViewResults” caption=“View Results?” control=“select” default=“No”>

          <attributes options=“Yes,No”/>

          <callbacks>

            <onvalidate>onViewPIDResults</onvalidate>

          </callbacks>

        </property>

      </propertygroup>

    </load>

  </simdata>

</extension>

Describing this in detail would take far longer than I have time for now, so I’m going to direct you to the ACT documentation. The gist of it is fairly simple though. XML provides a structured, hierarchical mechanism for describing the layout of the UI. Nested structures appear as child widgets of their parents. Callbacks are used within ACT to provide the hooks into the UI events so that we can respond to various user interactions. Beyond that, read the docs!! And, hey, before I hear any whining remember that in the old days of Jscript customization there wasn’t any documentation! When I was a Workbench Customization Kid we had to walk uphill, both ways, barefoot, in 8’ of snow for 35 miles… So shut it!

Making the Magic Happen

So, the UI is snazzy and all, but the heavy lifting really happens under the hood. Ultimately, what ACT provides us, when we are creating new BCs for ANSYS, is a clever way to insert APDL commands into the ds.dat input stream. Remember that at its core all Mechanical is, is a glorified APDL generator. I’m sure the developers love me reducing their hard labor to such mundane terms, but it is what it is… So, at the end of the day, our little ACT load objects are nothing more than miniature APDL writers. You thought we were doing something cool…

So, the magic happens when we collect up all of the input data the user specifies in our snazzy UI and turn that into APDL code that implemented the PID controller. This is obviously why I started by developing the APDL code first outside of WB. The APDL code is the true magic. Collecting up the user inputs and writing them to the ds.dat file occurs inside the getcommands callback. If you look closely at the XML code, you will notice two getcommands callbacks. The first one calls a python function named: writePIDThermostatLoad. This callback is scheduled to fire when Mechanical is finished writing all of the standard loads and boundary conditions that it implements natively and is about to write the first APDL solve command. Our commands will end up in the ds.dat file right at this location. I chose this location for the following reason. Our APDL code for the PID thermostat will be generating new element types and new nodes and elements not generated by Workbench. Sometimes workbench implements certain boundary conditions using surface effect element types. So, these native loads and boundary conditions themselves may generate new elements and element types. Workbench knows about those, because it’s generating them directly; however, it has no idea what I might be up to in my PID thermostat load. So, if it were to write additional boundary conditions after my PID load, it very well might reuse some of my element type ids, and even node/element ids. The safer thing to do is to write our APDL code so that it is robust in the presence of an unknown maximum element type/real constant set/node number/etc… Then, we insert it right before the solve command, after WB has written all of its loads and boundary conditions. Thus, the likelihood of id collisions is greatly reduced or eliminated.

Note, too, that ACT provides some utility functions to generate a new element type id and increment the internal counter within Workbench; however, I have found that these functions do not account for loads and boundary conditions. Therefore, in my testing thus far, I haven’t found them safe to use.

The second getcommands callback is setup to fire when Workbench finishes writing all of the solve commands and has moved to the post processing part of the input stream. I chose to implement a graphing functionality for displaying the relevant output data from the PID elements. Thus, I needed to retrieve this data from ANSYS after the solution is complete so that I can present it to the user. I accomplished this by writing a little bit of APDL code to enter /post26 and export all of the data I wish to plot to a CSV file. By specifying this second getcommands callback, I could instruct Workbench to insert the APDL commands after the solve completed.

Viewing the Results

Once the solution has completed, clicking on the “View Results?” dropdown and choosing “Yes” will bring up the following result viewer I implemented for the PID controller. All of the graphing functionality is provided by ACT in an import module called “chart”. This result viewer is simply implemented as a dialog with a single child control that is the ACT chart widget. This widget also allows you to layout multiple charts in a grid, as we have here. As you can see, we can display all of the relevant output data for the controls cleanly and efficiently using ACT! While this can also be accomplished in ANSYS Mechanical APDL, I think we would all agree that the results are far more pleasing visually when implemented in ACT.

clip_image006

Where Do We Go from Here?

Now that I’ve written an ACT module, my next steps are to clean it up and try to make it a little more production ready. Once I’m satisfied with it, I’ll publish it on this blog and on the appropriate ANSYS library. Look for more posts along the way if I uncover additional insights or gotchas with ACT programming. I will leave you with this, however. If you have put off ACT programming you really should reconsider! Being mostly new to ACT, I was able to get this little boundary condition hooked up and functioning in less than a week’s time. Given the way the user interface turned out and the flexibility thus far of the control, I’m quite pleased with that. Without the documentation and general availability of ACT, this effort would have been far more intense. So, try out ACT! You won’t be disappointed.

ACT Extension for a PID Thermostat Controller (PART 1)

I’m going to embark on a multipart blog series chronicling my efforts in writing a PID Thermostat control boundary condition for workbench. I picked this boundary condition for a few of reasons:

  1. As far as I know, it doesn’t exist in WB proper.
  2. It involves some techniques and element types in ANSYS Mechanical APDL that are not immediately intuitive to most users. Namely, we will be using the Combin37 element type to manage the control.
  3. There are a number of different options and parameters that will be used to populate the boundary condition with data, and this affords an opportunity to explore many of the GUI items exposed in ACT.

This first posting goes over how to model a PID controller in ANSYS Mechanical APDL.  In future articles I will share my efforts to refine the model and us ACT to include it in ANSYS Workbench.

PID Controller Background

Let’s begin with a little background on PID controllers. Full disclaimer, I’m not controls engineer, so take this info for what it is worth. PID stands for Proportional Integral Differential controller. The idea is fairly simple. Assume you have some output quantity you wish to control by varying some input value. That is, you have a known curve in time that represents what you would like the output to look like. For example:

image

The trick is to figure out what the input needs to look like in time so that you get the desired output. One way to do that is to use feedback. That is, you measure the current output value at some time, t, and you compare that to what the desired output should be at that time, t. If there is no difference in the measured value and the desired value, then you know whatever you have set the input to be, it is correct at least for this point in time. So, maybe it will be correct for the next moment in time. Let’s all hope…

However, chances are, there is some difference between what is measured and what is desired. For future reference we will call this the error term. The secret sauce is what to do with that information? To make things more concrete, we will ground our discussion in the thermal world and imagine we are trying to maintain something at a prescribed temperature. When the actual temperature of the device is lower than the desired temperature, we will define that as a positive error. Thus, I’m cold; I want to be warmer: that equals positive error. The converse is true. I’m hot; I want to be colder: that equals negative error.

One simple way we could try to control the input would be to say, “Let’s make the input proportional to the error term.” So, when the error term is positive, and thus I’m cold and wish to be warmer, we will add energy proportionate to the magnitude of the error term. Obviously the flip side is also true. If I’m hot and I wish to be cooler my negative error term would mean that remove energy proportionate to the magnitude of the error term. This sounds great! What more do you need? Well, what happens if I’m trying to hold a fixed temperature for a long time? If the system is not perfectly adiabatic, we still have to supply some energy to make up for whatever the system is losing to the surroundings. Obviously, this energy loss occurs even with the system is in a steady state configuration and at the prescribed temperature! But, if the system is exactly at the prescribed temperature, then the error term is zero. Anything proportionate to zero is… zero. That’s a bummer. I need something that won’t go to zero when my error term goes to zero.

What if I could keep a record of what I’ve done in the past? What if I accumulated all of the past error from forever? Obviously, this has the chance of being nonzero even if instantaneously my error term is zero. This sounds promising. Integrating a function of time with respect to time is analogous to accumulating the function values from history past. Thus, what if I integrated my error term and then made my input also proportional to that value? Wouldn’t that help the steady state issue above? Sure it would. Unfortunately, it also means I might go racing right on by my set point and it might take a while for that “mistake” to wash out of the system. Nothing is free. So, now I have kept a record of my entire past and used that to help me in the present, what if I could read the future? What if could extrapolate out in time?

Derivatives allow us to make a local extrapolation (in either direction) about a curve at a fixed point. So, if my curve is a function of time, which in our case the curves are, forward extrapolation is basically jumping ahead into the future. However, we can’t truly predict the future, we can only extrapolate on what has recently happened and make the leap of faith that it will continue to happen just as it has. So, if I take the derivative of my error term with respect to time, I can roll the dice a little a make some of my input proportional to this derivative term. That is, I can extrapolate out in time. If I do it right, I can actually make the system settle out a little faster. Remember that when the error term goes to zero and stays there, the derivative of the error term also goes to zero. So, when we are right on top of our prescribed value this term has no bearing on our input.

So, a PID controller simply takes the three concepts of how to specify an input value based on a function of the error term and mixes them together with differing proportions to arrive at the new value for the input. By “tuning” the system we can make it such that it responds quickly to change and it doesn’t wildly overshoot or oscillate.

Implementing a PID controller in ANSYS MAPDL

We will begin by implementing a PID controller in MAPDL before moving on to implementing the boundary condition in ANSYS Workbench via the ACT. We would like the boundary condition to have the following features:

  1. Ultimately we would like to “connect” this boundary condition to any number of nodes in our model. That is, we may want to have our energy input occur on a vertex, edge or face of the model in Workbench. So, we would like the boundary condition to support connecting to any number of nodes in the model.
  2. Likewise, we would like the “measured output” to be influenced by any number of nodes in our model. That is, if more than one node is included in the “measured value” set, we would like ANSYS to use the average temperature of the nodes in that set as our “measured output”. Again, this will allow us to specify a vertex, edge, face or body of the model to function as our measurement location. The measured value should be the average temperature on this entity. Averaging needs to be intelligent. We need to weight the average based on some measure that accounts for the relative influence of a node attached to large elements vs one attached to small elements.
  3. We would like to be able to independently control the proportional, integral and derivative components of the control algorithm.
  4. It would be nice to be able to specify whether this boundary condition can only add energy, only remove energy or if it can do both.
  5. We would like to allow the set point value to also be a function of time so that it too can change with time.
  6. Finally, it would be nice to be able to post process some of the heat flow quantities, temperature values, etc… associated with this boundary condition.

This is a pretty exhaustive list of requirements for the boundary condition. Fortunately, ANSYS MAPDL has built into it an element type that is perfectly suited for this type of control. That element type is the combin37.

Introducing the Combin37 Element Type

Understanding the combin37 element in ANSYS MAPDL takes a bit of a Zen state of mind… It’s, well, an element only a mother could love. Here is a picture lifted from the help:

clip_image003

OK. Clear as mud right? Actually, this thing can act as a thermostat whether you believe me from the picture or not. Like most/all ANSYS elements that can function in multiple roles, the combin37 is expressed in its structural configuration. It is up to you and me to mentally map it to a different set of physics. So, just trust me that you can forget the damping and FSLIDE and little springy looking thing in the picture. All we are going to worry about is the AFORCE thing. Mentally replace AFORCE with heat flow.

Notice those two little nodes hanging out there all by their lonesome selves labeled “control nodes”. I think they should have joysticks beside them personally, but ANSYS didn’t ask me. Those little guys are appropriately named. One of them, NODE K actually, will function as our set point node. That is, whatever temperature value we specify in time for NODE K, that same value represents the set point temperature we would like our “measured” location take on in time as well. So, that means we need to drive NODE K with our set point curve. That should be easy enough. Just apply a temperature boundary condition that is a function of time to that node and we’re good to go. Likewise, NODE L represents the “measured” temperature somewhere else in the model. So, we need to somehow hook NODE L up to our set of measurement nodes so that it magically takes on the average value of those nodes. More on that trick later.

Now, internally the combin37 subtracts the temperature at NODE K from NODE L to obtain an error term. Moreover, it allows us to specify different mathematical operations we can perform on the error term, and it allows us to take the output from those mathematical operations and drive the magical AFORCE thingy, which is heat flow. Guess what those mathematical operations are? If you guessed simply making the heat flow through the element proportional to the error, proportional to the time integral of the error and proportional to the time derivative of the error you would be right. Sounds like a PID controller doesn’t it? Now, the hard part is making sense of all the options and hooking it all up correctly. Let’s focus on the options first.

Key Option One and the Magic Control Value

Key option 1 for the combin37 controls what mathematical operation we are going to perform on the error term. In order to implement a full PID controller, we are going to need three combin37 elements in our model with each one keyed to a different mathematical operation. ANSYS calls the result of the mathematical operation, Cpar. So, we have the following:

KEYOPT(1) Value Mathematical Operation
0,1 image
2 image
3 image
4 image
5 image

Thus, for our purposes, we need to set keyopt(1) equal to 1,4 and 2 for each of the three elements respectively.

Feedback is realized by taking the control parameter Cpar and using it to modify the heat flow through the element, which is called AFORCE. The AFORCE value is specified as a real constant for the element; however, you can also rig up the element so that the value of AFORCE changes with respect to the control parameter. You do this by setting keyopt(6)=6. The manner in which ANSYS adjusts the AFORCE value, which again is heat flow, is described by the following equation:

image

Thus, the proportionality constant for the Proportional, Integral and Derivative components are specified with the C1 variable. RCONST, C3 and C4 are all set to zero. C2 is set to 1. Also note that ANSYS first takes the absolute value of the control parameter Cpar before plugging it into this equation. Furthermore, the direction of the AFORCE component is important. A positive value for AFORCE means that the element generates an element force (heatflow) in the direction specified in the diagram.  That is, it acts as a heat sink. So, assuming the model is attached to node J, the element acts as a heat sink when AFORCE is positive. Conversely, when AFORCE is negative, the element acts like a heat source. However, due to the absolute value, Cpar can never take on a negative value. Thus, when this element needs to act as an energy source to add heat to our model, the coefficient C1 must be negative. The opposite is true when the element needs to act as an energy sink.

Key Option Four and Five and when is it Alive?

If things weren’t confusing enough already, hold on as we discuss Keyopt 4 and 5. Consider the figure below, again lifted straight from the help.

clip_image016

The combination of these two key options controls when the element switches on and becomes “alive”. Let’s take the simple case first. Let’s assume that we are adding energy to the model in order to bring it up to a given temperature. In this case, Cpar will be positive because our set point is higher than our current value. If the element is functioning as a heat source we would like it to be on in this condition. Furthermore, we would like it to stay on as long as our error is positive so that we continue adding energy to bring the system up to temperature. Consider the diagram in the upper left. Imagine that we set ONVAL = 0 and OFFVAL = 0.001. Whenever Cpar is greater than ONVAL.  So this sounds like exactly what we want when the element is functioning as a heat source. Thus, keyopt(4)=0 and keyopt(5)=0.001 with OFFVAL=ONVAL=0 is what we want when the element needs to function as a heat source.

What about when it is a heat sink?  In this case we want the element to be active when the error term is negative; that is, when the current temperature is higher than the set point temperature.  Consider the diagram in the middle left.  This time let OFFVAL=0 and OFFVAL=-0.001.  In this case, whenever Cpar is negative (less than OFFVAL) then the element will be active.  Thus, keyopt(4)=0 and keyopt(5)=1 with OFFVAL=-0.001 ONVAL=0 is what we want when the element needs to function as a heat sink.  Notice, that if you set ONVAL=OFFVAL then the element will always stay on; thus, we need to provide the small window to activate the switching nature of the element.

Thus, we see that we need six different combin37 elements, three for a PID controlled heat sink and three for a PID controlled heat source, to fully specify a PID controlled thermal boundary condition. Phew… However, if we set all of the proportionality constants for either set of elements defining the heat sink or heat source to zero, we can effectively turn the boundary condition into only a heat source or only a heat sink, thus meeting requirement four listed above. While we’re marking off requirements, we can also mark off requirements three and five. That is, with this combination of elements we can independently control the P, I and D proportionality constants for the controller. Likewise, by putting a time varying temperature constraint on control node K, we can effectively cause the set point value to change in time. Let’s see if we can now address requirements one and two.

How do we Hook the Combin37 to the Rest of the Model?

We will address this question in two parts. First, how do we hook the “business” end of the combin37 to the part of the model to which we are going to add or remove energy? Second, how do we hook the “control” end of the combin37 to the nodes we want to monitor?

Hooking to the Combin37 to the Nodes that Add or Remove Energy

To hook the combin37 to the model so that we can add or remove energy we will use the convection link elements, link34. These elements effectively act like little thermal resistors with the resistance equation being specified as:

image

In order to make things nice, we need to “match” the resistances so that each node effectively sees the same resistance back to the combin37 element. We do this by varying the “area” associated with each of these convective links. To get the area associated with a node we use the arnode() intrinsic function. See the listing for details.

Hooking the Combin37 to the Nodes that Function as the Measured Value

As we mentioned in our requirements, we would like to be able to specify more than one or more nodes to function as the measured control value for our boundary condition. More precisely, if more than one node is included in the measurement node set, we would like ANSYS to average the temperatures at those nodes and use that average value as the measurement temperature. This will allow us to specify, for example, the average temperature of a body as the measurement value, not just one node on the body somewhere. However, we would also like for the scheme to work seamlessly if only one node is specified. So, how can we accomplish this? Constraint equations come to our rescue.

Remember that a constraint equation is defined as:

image

How can we use this to compute the average temperature of a set of nodes, and tie the control node of the combin37 to this average? Let’s begin by formulating an equation for the average temperature of a set of nodes. We would like this average to not be simply a uniform average, but rather be weighted by the relative contribution a given node should play in the overall average of a geometric entity. For example, assume we are interested in calculating the average temperature of a surface in our model. Obviously this surface will have associated with it many nodes connected to many different elements. Assume for the moment that we are interested in one node on this face that is connected to many large elements that span most of the area of this face. Shouldn’t this node’s temperature have a larger contribution to the “average” temperature of the face as say a node connected to a few tiny elements? If we just add up the temperature values and divide by the number of nodes, each node’s temperature has equal weight in the average. A better solution would be to area weight the nodal temperatures based on the area associated with each individual node. Something like:

image

That looks a little like our constraint equation. However, in the constraint equation I have to specify the constant term, whereas in the equation above, that is the value (Tavg) that I am interested in computing. What can I do? Well, let’s add in another node to our constraint equation that represents the “answer”. For convenience, we’ll make this the control node on our combin37 elements since we need the average temperature of the face to be driving that node anyway. Consider:

image

Now, our constant term is zero, and our Ci’s are Ai/AT and -1 for the control node. Voila! With this one constraint equation we’ve compute an area weighted average of the temperature over a set of nodes and assigned that value to our control node. CE’s rock!

An Example Model

This post is already way too long, so let’s wrap things up with a little example model. This model will illustrate a simple PI heat source attached to an edge of a plate with a hole. The other outer edges of the plate are given a convective boundary condition to simulate removing heat. The initial condition of the plate is set to 20C. The set point for the thermostat is set to 100C. No attempt is made to tune the PI controller in this example, so you can clearly see the effects of the overshoot due to the integral component being large. However, you can also see how the average temperature eventually settles down to exactly the set point value. clip_image026

The red squiggly represents where heat is being added with the PI controller. The blue squiggly represents where heat is being removed due to convection. Here is a plot of the average temperature of the body with respect to time where you can see the response of the system to the PI control.

clip_image027

Here is another run, where the set point value ramps up as well. I’ve also tweaked the control values a little to mitigate some of the overshoot. This is looking kind of promising, and it is fun to play with. Next time we will look to integrate it into the workbench environment via an actual ACT extension.

clip_image028

Part 2 is here

Model Listing

I’ve included the model listing below so that you can play with this yourself. In future posts, I will elaborate more on this technique and also look to integrate it into an ACT module.

 

finish

/clear

 

/prep7

*get,etmax_,etyp,0,num,max

P_et=etmax_+1

I_et=etmax_+2

D_et=etmax_+3

Link_et=etmax_+4

mass_et=etmax_+5

 

 

 

et,P_et,combin37

et,I_et,combin37

et,D_et,combin37

et,Link_et,link34

et,mass_et,mass71

Kp=1

Ki=2

Kd=0

 

keyopt,P_et,1,0    ! Control on UK-UL

keyopt,P_et,2,8    ! Control node DOF is Temp

keyopt,P_et,3,8    ! Active node DOF is Temp

keyopt,P_et,4,0    ! Wierdness for the ON/OFF range

keyopt,P_et,5,0    ! More wierdness for the ON/OFF range

keyopt,P_et,6,6    ! Use the force, Luke (aka Combin37)

keyopt,P_et,9,0    ! Use the equation, Duke (where is Daisy…)

 

 

keyopt,I_et,1,4    ! Control on integral wrt time

keyopt,I_et,2,8    ! Control node DOF is Temp

keyopt,I_et,3,8    ! Active node DOF is Temp

keyopt,I_et,4,0    ! Wierdness for the ON/OFF range

keyopt,I_et,5,0    ! More wierdness for the ON/OFF range

keyopt,I_et,6,6    ! Use the force, Luke (aka Combin37)

keyopt,I_et,9,0    ! Use the equation, Duke (where is Daisy…)

 

 

keyopt,D_et,1,2    ! Control on first derivative wrt time

keyopt,D_et,2,8    ! Control node DOF is Temp

keyopt,D_et,3,8    ! Active node DOF is Temp

keyopt,D_et,4,0    ! Wierdness for the ON/OFF range

keyopt,D_et,5,0    ! More wierdness for the ON/OFF range

keyopt,D_et,6,6    ! Use the force, Luke (aka Combin37)

keyopt,D_et,9,0    ! Use the equation, Duke (where is Daisy…)

 

keyopt,mass_et,3,1 ! Interpret real constant as DENS*C*Volume

 

 

 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!        S M A L L   T E S T   M O D E L       !!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

test_et=etmax_+10

et,test_et,plane55

 

mp,kxx,test_et,70

mp,dens,test_et,8050

mp,c,test_et,0.4

! Thickness of plate

r,test_et,0.1

 

! Plane55 element

keyopt,test_et,3,3

 

! Make a block

k,1

k,2,1,0

k,3,1,1

k,4,0,1

a,1,2,3,4

! Make a hole

cyl4,0.5,0.5,0.25

! Punch a hole

asba,1,2

 

type,test_et

mat,test_et

real,test_et

esize,0.025

amesh,all

 

! create an nodal component for the

! ‘attachment’ location

lsel,s,loc,x,0

nsll,s,1

cm,pid_attach_n,node

 

! create a nodal component for the

! ‘monitor’ location

allsel,all

cm,pid_monitor_n,node

 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!        B E G I N   P I D   M O D E L         !!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 

! Real constant and mat prop for the mass element

mp,qrate,mass_et,0 ! Zero heat generation rate for the element

r,mass_et,1e-10  ! Extremely small thermal capacitance

 

! Material properties for convection element

! make the convection “large”

mp,hf,Link_et,1e5

 

! Real constant for the combin37 elements

! that ack as heaters

on_val_=0

off_val_=1e-3

r,P_et,0,0,0,on_val_,off_val_,0

rmore,0,1,-Kp,1

 

r,I_et,0,0,0,on_val_,off_val_,0

rmore,0,1,-Ki,1

 

r,D_et,0,0,0,on_val_,off_val_,0

rmore,0,1,-Kd,1

 

! build the PID elements

*get,nmax_,node,0,num,max

BASE_NODE=nmax_+1

P_NODE_J=nmax_+2

I_NODE_J=nmax_+3

D_NODE_J=nmax_+4

NODE_K=nmax_+5

NODE_L=nmax_+6

! Create the nodes.  They can be all coincident

! as we will refer to them solely by their number.

! They will be located at the origin

*do,i_,1,6

    n,nmax_+i_

*enddo

 

! Put a thermal mass on the K and L nodes

! for each control element to give them

! thermal DOFs

type,mass_et

mat,mass_et

real,mass_et

e,NODE_K

e,NODE_L

 

! Proportional element

type,P_et

mat,P_et

real,P_et

e,BASE_NODE,P_NODE_J,NODE_K,NODE_L

 

! Integral element

type,I_et

mat,I_et

real,I_et

e,BASE_NODE,I_NODE_J,NODE_K,NODE_L

 

! Derivative Element

type,D_et

mat,D_et

real,D_et

e,BASE_NODE,D_NODE_J,NODE_K,NODE_L

 

 

 

! Ground the base node

d,BASE_NODE,temp,0

 

! Get a list of the attachment nodes

cmsel,s,pid_attach_n

*get,numnod_,node,0,count

attlist_=

*dim,attlist_,,numnod_

*vget,attlist_(1),node,0,nlist

*get,rmax_,rcon,0,num,max

 

! Hook the attachment nodes to the

! end of the control element with

! convection links

*do,i_,1,numnod_

    n1_=attlist_(i_)

    a1_=arnode(n1_)

    r1_=rmax_+i_

    r,r1_,a1_

    type,link_et

    mat,link_et

    real,r1_

    e,P_NODE_J,n1_

    e,I_NODE_J,n1_

    e,D_NODE_J,n1_

*enddo

 

! Hook up the monitor nodes

cmsel,s,pid_monitor_n

*get,numnod_,node,0,count

monlist_=

monarea_=

*dim,monlist_,,numnod_

*dim,monarea_,,numnod_

*vget,monlist_(1),node,0,nlist

! We are going to need these areas

! so, hold on to them

*do,i_,1,numnod_

    monarea_(i_)=arnode(i_)

*enddo

*vscfun,totarea_,sum,monarea_(1)

! Write the constraint equations

/pbc,ce,,0

ce,next,0,NODE_L,temp,-1

*do,i_,1,numnod_

    ce,high,0,monlist_(i_),temp,monarea_(i_)/totarea_

*enddo

 

! Create a transient setpoint temperature

*dim,setpoint_,table,2,,,time

setpoint_(1,0)=0

setpoint_(2,0)=3600

setpoint_(1,1)=100

setpoint_(2,1)=100

 

! Constrain the temperature node to be

! at the setpoint value

allsel,all

d,NODE_K,temp,%setpoint_%

 

 

! Apply an initial condition of

! 20 C to everything

allsel,all

ic,all,temp,20

 

/solu

antype,trans,new

time,1000

deltim,1,1,1

lsel,s,loc,x,0.1,1

nsll,s,1

sf,all,conv,20,20

 

allsel,all

outres,all,all

solve

/post26

! Plot the response temperature

! and the setpoint temperature

nsol,2,NODE_L,temp,,temp_r

nsol,3,NODE_K,temp,,temp_s

xvar,1

plvar,2,3

ANSYS Mechanical Scripting: HOWTO Part 5

Where we Left Off Last Time…

In the last post of this series, I left you with an introduction to the Microsoft Visual Studio debugger and a description of how to debug scripts that we write.  In this post, I’m going to illustrate how to write what might actually turn out to be a useful script!  The idea for this script came from a request from one of our users, Craig Hildreth, who has always wished that there was a better way to define slice planes in Mechanical. 

Problem Description

Lets say you want to create a slice plane in ANSYS Mechanical that cuts precisely through a certain location in your model and that is precisely oriented to be normal to a certain direction.  How would you do it?  If you’ve used the slice tool in Mechanical before, you know that you basically yank a line segment with the mouse across the screen.  The resulting slice plane is defined such that if you were to extend the line segment you created with the mouse to infinity and then extrude the line in the direction of the current view the resulting plane formed by this operation would be the slice plane.  Can it be more intuitive?  So, of you go spinning and panning your model around to where you think you are looking right down the slice plane, then with the steadiness of a surgeon you try with all your might to draw a straight line with the mouse.  (Good luck!)  Now, to ANSYS’ credit, this is a really cool tool if you just want to whip a quick slice plane on the model to look inside somewhere.  Especially when you’ve got that cool little drag widget.  However, if you are trying to use the slice tool for any kind of meaningful post processing, your probably out of luck.  What would be awfully darn useful would be create or use a Cartesian coordinate system that would then define the slice plane.  That sounds kind of intuitive, right?  Anyway, you might optimistically think that functionality should already be there, but alas to my knowledge it is not.  So lets see if we can fake it.  (If you just want the code and don’t care how we write it, check out this post where we give the code.  This is a long post because it turns out this is trickier than it looks…)

How the Heck Do You Create a Slice Plane… Programmatically that is?

So, step one is to figure out what Mechanical does under the hood when we create a slice plane.  If you haven’t read Part 2of this series, I suggest you do so now.  We’re going to use the techniques in that post to figure out what Mechanical is doing when we create a slice plane. 

Searching the ANSYS Source Files for Info on Creating a Slice Plane

If you haven’t started up Mechanical, do so now and make sure that the slice plane window is visible by ensuring that menu item View->Windows->Section Planes has a check mark beside it.  In the section plane window, hover your mouse of the create section plane button.  You should see a tool tip like the following:

image

We see that the tooltip is called “New Section Plane”.  This is a string for which we can search.  So, search the entire directory tree C:\Program Files\ANSYS Inc\v130\aisol\DesignSpace\ for the string “New Section Plane”.  We see that in C:\Program Files\ANSYS Inc\v130\aisol\DesignSpace\DSPages\Language\en-us\xml\dsstringtable.xml we find that string.  Remember that this is where you are normally going to find GUI related strings.  This string has the string id: ID_NewSectionPlane.  So, lets search for that.  One place we find this string is inside the file: C:\Program Files\ANSYS Inc\v130\aisol\DesignSpace\DSPages\scripts\DSSectionPlanesScript.js  This file name sounds promising.  Here is some code from that file showing where the id shows up.

 

btnAdd = toolbar.Buttons.AddButton( TBT_BUTTON, 
"NewSectionPlane", ID_ADD, 
"NewSectionPlane", false, "", 
localString("ID_NewSectionPlane") );

btnAdd.OnClick.AddCOM(g_sectionPaneObjects, 
"SectionPane_Add");
  

This code is reformatted to fit better in the blog.  You can see that the code appears to be adding a button to a toolbar and also registering a callback for when the button is clicked.  The callback is what we’re interested in, so lets search for “SectionPane_Add”. Searching for that leads us to the following code:

this.SectionPane_Add = function sectionPane_Add(sender, args)
{
    if (!AllowSlicing ())  return;
    doAddSectionPlane();
 }

We’re getting closer!  It looks like we need to peek into doAddSectionPlane().  Let’s search for that.  It turns out that that function is a bit of a dead end, but lucky for us right below it we see the following code:

function doDeleteSectionPlane()
{
    var sliceTool  = ds.Graphics.SliceTool;
    if (sliceTool == null)
        return;
   ...
}

Hey!  There is something called a SliceTool that is a member of the Graphics object associated with the global DesignSpace Object.  Lets quit grepping and slogging through the ANSYS Mechanical code base and fire up the debugger to see what we can learn about the SliceTool.

Using the Visual Studio Debugger to Learn About the Slice Tool

Now that we’ve figured out that there is an object called SliceTool that is a member of the Graphics object, lets use the debugger to interrogate that object.  Create the following javascript code and save it somewhere on your hard drive.

function slice_tool_interrogate() {
    debugger;
    var slice_tool = DS.Graphics.SliceTool;
    // Useless code to give us a place
    // to stop
    var i = 1;
}

slice_tool_interrogate();

The code above will allow us to stop in the debugger and look at the slice tool object.  The following screen shot shows the slice_tool object’s properties and methods in the debugger.

image

Hmm… Well, there is one function there called SetPlane(ix1, iy1, ix2, iy2) that looks kind of promising.  Unfortunately, I don’t see any function that would seem to hint at defining the plane with a point and normal, for example.  It looks like the programmatic interface mimics the user interface.  Arg…  So, lets just guess that given two points in screen coordinates (x1, y1 and x2,y2) we can create a slice plane programmatically.  What should x1, y1 and x2, y2 be?  Are screen coordinates normalized? Where is the screen origin; top left or bottom right or middle of the screen? The fact that the arguments have an “I” in front of them hints that perhaps they are integers, which would rule out normalized values.  What to do?  Well, lets try this.  Lets orient some geometry so that it is mostly centered with the screen and, using the immediate window in the debugger, pass in some values to see what happens.  Here is a screen shot prior to using the immediate window to test things out.

image

Now, I’m going to run our little test script again under the debugger, but after the slice_tool object is created, I’m going to enter the following into the debugger’s immediate window:

slice_tool.SetPlane(0,0.5,1,0.5)

When I do that, I get the following:

image

So, we created a slice plane!!!  Unfortunately, it doesn’t appear that it cut through anything.  We learned a couple of things.

  1. The function DS.Graphics.SliceTool.SetPlane(x1, y1, x2, y2) will indeed create a slice plane object.
  2. The coordinates don’t appear to be normalized.  (Can’t rule this out definitively, but lets try something else.)

Delete this plane, and lets work with the hypothesis that we need to pass in the screen coordinates as integer values.  If I were a programmer, I would choose pixels as my unit of measure for screen coordinates.  Things like mouse clicks in particular can be directly translated easily.  So, now the question is how big is the screen?  Or more precisely, how big is the graphics window?  This DS.Graphics object might hold the answers.  If you poke around this object inside the debugger you’ll notice that it has a property called GfxWindow.  That sounds promising!  If we look inside that object, we see:

image 

Ah ha!  This thing has a height and a width property associated with it, which I assume by the values associated with these properties are in units of pixels.  So lets try the following code.  You can modify your little test script as follows:

function slice_tool_interrogate() {
    debugger;
    var slice_tool = DS.Graphics.SliceTool;
    var height = DS.Graphics.GfxWindow.height;
    var width = DS.Graphics.GfxWindow.width;
    var slice_plane = 
        slice_tool.SetPlane(0, height / 2, width, height / 2);
    // Useless code to give us a place
    // to stop
    var i = 1;
}
slice_tool_interrogate();

When I run this little script.  My “after photo” looks like:

image

Ha!  We got ourselves a genuine slice plane.  And, we got it where we want it.  Well, sort of…  OK, so forgive me for celebrating a minor victory, but we were able to create a slice plane and draw a perfectly straight line right across the middle of the screen.   Now, if only we could programmatically rotate and pan the model around so that we are centered right on the origin of the plane we are interested in and we are looking right down parallel to the plane.  If we can do that, then we just draw a quick line with the SetPlane function and voila, we’re done.

How the Heck to You Rotate and Pan the View… Programmatically that is?

You know, this Graphics object is pretty handy.  Lets see what other goodies it might hold.  If we poke around in it a bit, we see that it has a property called Camera.  Camera sounds like something that would help us look at things from a certain point of view.   Here is a peek into the camera object inside the debugger:

image

This is looking promising.  You’ll notice that we’ve got a pan and rotate function.  And, it looks like we’ve got a focus point and view direction vector.  This is looking good.  Maybe we can write a function that will pan the view around until the origin of our coordinate system is centered where we want it and simultaneously rotate the view around until our view direction is oriented how we want it.  One other thing we’ll have to figure out is that we could rotate the view about the view direction vector so that it spins around the axis we are looking down.  Think of a drill bit spinning.  In that case neither the view direction or the focus point would likely change, but there is that rotateZ() function that might bail us out.  Lets see if we can figure out what all these things do/mean.  Lets start with the following test functions:

function do_pan(up, right) {
    var camera = DS.Graphics.Camera;
    var graphics = DS.Graphics;
    var a = 0;
    for (var i = 0; i < 50; i++) {
        camera.pan(up, right);
        // Just spin her so that the "animation" is more smooth.
        for (var j = 0; j < 1000; j++) {
            a += Math.sin(j);
        }
        graphics.Refresh();
    }
    // I don't know if Javascript does dead code removal, but
    // try to trick any optimizations the javascript interpreter
    // might do to our busy loop.
    return a;
}

function do_rotate(up, right) {
    var camera = DS.Graphics.Camera;
    var graphics = DS.Graphics;
    var a = 0;
    for (var i = 0; i < 50; i++) {
        camera.rotate(up, right);
        // Just spin her so that the "animation" is more smooth.
        for (var j = 0; j < 1000; j++) {
            a += Math.sin(j);
        }
        graphics.Refresh();
    }
    // I don't know if Javascript does dead code removal, but
    // try to trick any optimizations the javascript interpreter
    // might do to our busy loop.
    return a;
}

These two functions basically call the Camera pan and rotate functions repeatedly with some values passed in.  Note that if you include these in your little test script, once you stop in the debugger, you can call these functions inside the immediate window and watch what happens.  Just type do_pan(5,0) for example.  If you play around with these functions, you’ll realize that the rotate function seems to rotate the model about either the screen vertical axis or the screen horizontal axis.  Interestingly, the pan function seems to be backwards.  That is, the amountUp argument seems to move the model to the right and the amountRight argument seems to move the model up.  Ah well…  We’ll just make note of it.  How do we make sense of all these things?

Math, Math… Wonderful Math! Vectors and Rotations Galore!

Note, if math makes your head spin, we’re going to spin.  We’re going to take dot products and cross products and do all kinds of other gyrations.  So, feel free to skip this section if you are not interested.  However, this section constitutes the meat of figuring out what Mechanical is doing with the graphics given the few bits of information we have highlighted above. 

Lets start with rotating the view so that we’re looking where we want to look.  Our rotate function spins the model about the screen vertical and horizontal axes, so I’m going to propose an iterative technique that will sequentially rotate up and rotate right until the view vector is parallel to a supplied orientation vector.  Remember that we need to look edge on to the plane we wish to use as a cut plane so that we can draw a straight line across our screen and create the cut plane.  So, our view vector needs to be “in the plane” of our cut plane.  Consider the image below:

view_direction

The red vector represents the current view direction vector.  Now, assume that we call the rotate function associated with the camera object and request an up rotation of 15 degrees.  The resulting new view direction vector is represented by the green vector above.  Lets call this the test vector.  Taking the cross product of these two vectors (green and red) allows us to determine the normal vector, shown in black above, that represents an axis about which our requested rotation occurred.  We now know that if we call the rotate function with a value for the up rotation given the current view vector, we will rotate about the normal vector above.  So, now rotate the view back to the original view direction.  Let’s assume that our desired view direction is show above as the blue vector.  In general, it will not be in the plane of our current rotation.  However, if we project this vector onto the plane defined by the normal vector and the origin we obtain the purple vector above.  Now I know that I can rotate the view in the “up” direction by some angle and get to the purple vector. Its not the blue vector that I want, but it is the projection of the view vector on the plane in which my rotations are occurring. I can calculate a good guess for this angle by taking the ratio of the angle between the green and red vector with respect to my original trial rotation of 15 degrees.  Then, if I take the angle between the purple vector and the red vector, divide by my ratio, I’ll get a good guess as to what I should try to rotate my view so that I get to the purple vector.  If it works according to plan, this will get me closer to my desired view direction.  Now, I repeat this whole process by rotating in the other direction, that is rotating to the right.  Between each of these update rotations, I check to see if my current view vector is parallel to my desired view vector within some small tolerance.  If things go well, I should walk my way around to the desired view vector.  Consider the code below.

// Calculate the dot product of two vectors
function dot(va, vb) {
    return va[0] * vb[0] + va[1] * vb[1] + va[2] * vb[2];
}

// Calculate the length of a vector
function norm(va) {
    return Math.sqrt(va[0] * va[0] + va[1] * va[1] + va[2] * va[2]);
}

// Normalize a vector
function normalize(va) {
    var len = norm(va);
    if (len != 0.0) {
        return [va[0] / len, va[1] / len, va[2] / len];
    }
}

// Calculate the cross product of two vectors
function cross(va, vb) {
    return [va[1] * vb[2] - va[2] * vb[1],
            va[2] * vb[0] - va[0] * vb[2], 
            va[0] * vb[1] - va[1] * vb[0]];
}

// Scale a vector by a scalar
function scale(va, sc) {
    return [va[0] * sc, va[1] * sc, va[2] * sc];
}

// Add two vectors
function add(va, vb) {
    return [va[0] + vb[0],
            va[1] + vb[1],
            va[2] + vb[2]];
}

// Subtract vector b from vector a
function sub(va, vb) {
    return [va[0] - vb[0],
            va[1] - vb[1],
            va[2] - vb[2]];
}

function are_parallel(va, vb, tol) {
    return 1.0 - (Math.abs(dot(va, vb)) / (norm(va) * norm(vb))) < tol ? true : false;
}

function are_perpendicular(va, vb, tol) {
    return (Math.abs(dot(va, vb)) / (norm(va) * norm(vb))) < tol ? true : false;
}

// Rotate the view around so that we are looking in the 
// direction of the desired view
function rotate_view(desired_view) {
    // Get the camera and graphics objects
    var camera = DS.Graphics.Camera;
    var graphics = DS.Graphics;
    // This is our tolerance for being parallel to the desired view
    var eps = 1e-7;
    // This is the maximum number of iterations we'll try.  
    // While loops with a tolerance check only scare me...
    var max_iter = 200;
    // Get the current view
    var view_c = normalize([camera.ViewDirectionX,
                            camera.ViewDirectionY,
                            camera.ViewDirectionZ]);
    // Make sure we're normalized
    dvd = normalize(desired_view);
    // This should be close to 1 if parallel.

    var cnt = 1;
    var b_first_arg = true;
    var normal = null;
    var trial = null;
    var view_p = null;
    var trial_ang = 15;
    var applied_rotation = null;
    var previous_up = 0;
    var previous_right = 0;
    var right_factor = up_factor = -1.0;
    // Loop until we're parallel, or we give up trying
    do {
        var factor = cnt / max_iter;
        factor *= factor;
        factor *= factor;
        // Get the view direction
        vd = normalize([camera.ViewDirectionX,
                        camera.ViewDirectionY,
                        camera.ViewDirectionZ]);
        // Get the current up vector
        camera.rotate(trial_ang, 0);
        trial = normalize([camera.ViewDirectionX,
                       camera.ViewDirectionY,
                       camera.ViewDirectionZ]);
        current_up = normalize(cross(trial, vd));
        // Rotate back so we don't lose our place
        camera.rotate(-trial_ang, 0);
        

        vd = normalize([camera.ViewDirectionX,
                        camera.ViewDirectionY,
                        camera.ViewDirectionZ]);

        if (b_first_arg) {
            camera.rotate(trial_ang, 0);
        } else {
            camera.rotate(0, trial_ang);
        }
        trial = normalize([camera.ViewDirectionX,
                           camera.ViewDirectionY,
                           camera.ViewDirectionZ]);
        n = normalize(cross(vd, trial));
        // Rotate back
        if (b_first_arg) {
            camera.rotate(-trial_ang, 0);
        } else {
            camera.rotate(0, -trial_ang);
        }
        vd = normalize([camera.ViewDirectionX,
                        camera.ViewDirectionY,
                        camera.ViewDirectionZ]);
        // Only do this rotation if our desired view vector
        // is not nearly parallel to our current axis of rotation
        // If it is, then the code inside the if statement will
        // be skipped.
        if (!are_parallel(dvd, n, eps)) {
            // Project the desired view vector onto our plane of rotation
            vp = normalize(cross(n, normalize(cross(dvd, n))));
            // Calculate the angle between the projected vector
            // and our current view direction
            dot_product = dot(vp, vd);
            needed_rotation = Math.acos(dot(vp, vd)) * 180 / Math.PI
            
            // Knock it down by a factor associated with our iteration
            // scheme.  This helps prevent spurious jittering as we 
            // make our way there.
            needed_rotation *= (1.0 - factor);
            if (b_first_arg) {
                // If we start to diverge, try rotating back the other way
                if(previous_up < needed_rotation) {
                    up_factor = (up_factor < 0.0) ? 1.0 : -1.0;
                    previous_up = needed_rotation;
                }
                camera.rotate(up_factor * needed_rotation, 0.0);
                
                
            } else {
                if (previous_right < needed_rotation) {
                    right_factor = (right_factor < 0.0) ? 1.0 : -1.0;
                    previous_right = needed_rotation;
                }
                camera.rotate(0.0, right_factor * needed_rotation);
            }
        }
        // See if we are there yet
        vd = normalize([camera.ViewDirectionX,
                        camera.ViewDirectionY,
                        camera.ViewDirectionZ]);
        b_rotate_converged = are_parallel(vd, dvd, eps);
        if (b_rotate_converged) break;
        // Flip directions
        b_first_arg = !b_first_arg;
//        DS.Graphics.Refresh();
    } while (cnt++ < max_iter);
        
}

You can try this function out with some test vectors to see if indeed it rotates the view into the desired orientation.  One thing you will notice, however, is that the view sometimes is spun at a funny angle as we look down this axis.  So all we’ve done is ensured that the camera is pointing in a certain direction, but there still is one degree of freedom, so to speak, left unspecified with the camera.  That is, you could spin the camera any way you want around this view direction and still technically be looking down the view direction.  So, now how do we figure out what “spin” has been applied to the camera?  Moreover, what spin do we want applied?

The answer for the second question above is related to the coordinate system that will define our cut plane. Lets assume that we are wanting to look down the Y axis of the coordinate system which will define our cutting plane, and we want the Z axis of this coordinate system to be pointing straight up on the screen.  If we do this, when we draw a straight cut line horizontally across the middle of the screen, we will be defining a cut plane in the XY plane of the specified coordinate system.  That is exactly what we want.  In the code above we’ve got the part that will enable us to look down the Y axis, now we just need to figure out how to make the Z axis point straight up

In most computer graphics implementations, there is this notion of an “Up Vector” which is defined to be perpendicular to the line of sight and is used to determine “which way is up”.  That is, the up vector points in the positive Y axis on the screen.  If you look back up at the screen shot for the camera object displayed in the debugger, you will indeed see that there is an UpVector as a member of that object!  Yippee!!! Unfortunately, it says “Bad Variable Type”.  Bummer.  That’s code for “this variable exists but you can’t access it in the scripting language”.  We do have this variable, AngleRotationZ, that looks promising, but what is the reference for measuring this angle?  Also, we have a function called rotateZ(), which again looks promising, but how are we going to know how much to rotate?  We need to know which way is up, and how much to rotate around the camera’s Z axis so that we can orient our view such that when we draw a straight line across the middle of the screen, we’re cutting the model exactly on our desired plane. 

Remember that the first argument to the function camera.rotate(amountUp, amountRight) will rotate the view about the Up Vector.  Consider the image below:

rotate_z

Lets assume that we start with our current view direction vector, which is represented in the image above by the red vector.  Next, we perform a trial rotation of 15 degrees about the unknown up vector using the camera.rotate(15,0).  This moves our view direction vector the green vector shown above.  By taking the cross product of these two vectors, we can determine which way is “up” in the world coordinate system.  We then rotate the view back to where we started.  Now that we know what the current up vector is, we can take the dot product with the desired up vector, shown in purple above, and determine the angle about the camera Z axis we need to rotate so that up is truly up.  Not that this assumes that the desired up vector is perpendicular to the current view direction.  Consider the code below.

function rotate_up(desired_up_vector) {
    
    // Get the camera and graphics objects
    var camera = DS.Graphics.Camera;
    var graphics = DS.Graphics;
    var view_c = null;
    var trial = null;
    var current_up = null;
    var theta = null;
    var trial_ang = 15;
    var theta1;
    var theta2;
    var eps = 0.01;
    // Make sure our passed in value is normalized
    duv = normalize(desired_up_vector);

    // Which way are we currently looking now?
    vd = normalize([camera.ViewDirectionX,
                    camera.ViewDirectionY,
                    camera.ViewDirectionZ]);

    // Perform our trial rotation about the screen vertical axis
    camera.rotate(trial_ang, 0);
    trial = normalize([camera.ViewDirectionX,
                       camera.ViewDirectionY,
                       camera.ViewDirectionZ]);
    up = normalize(cross(trial, vd));
    // Rotate back so we don't lose our place
    camera.rotate(-trial_ang, 0);
    // How much to we need to rotate?
    theta1 = Math.acos(dot(duv, up)) * 180 / Math.PI;
    // Rotate assuming a positive rotation
    camera.rotateZ(theta1);
    DS.Graphics.Refresh();
    // Now we need to test to see if this is indeed the proper direction
    // Which way are we currently looking now?
    vd = normalize([camera.ViewDirectionX,
                    camera.ViewDirectionY,
                    camera.ViewDirectionZ]);
    // Perform our trial rotation about the screen vertical axis
    camera.rotate(trial_ang, 0);
    trial = normalize([camera.ViewDirectionX,
                       camera.ViewDirectionY,
                       camera.ViewDirectionZ]);
    up = normalize(cross(trial, vd));
    // Rotate back so we don't lose our place
    camera.rotate(-trial_ang, 0);
    // How much to we need to rotate?
    theta2 = Math.acos(dot(duv, up)) * 180 / Math.PI;
    if (theta2 > theta1) {
        // Rotate backwards if our second guess is larger than our first
        camera.rotateZ(-theta2);
//        DS.Graphics.Refresh();
    }
}

In this code we implement the algorithm sketched out above.  The only little trick is that we might end up rotating in the wrong direction.  We can determine this by peforming the same test again, and if we are still off then we can rotate back the other direction.

Now that we can look in a particular direction and we can deduce which was is up, all that is left is to pan the view around so that the location of our cut plane is right in the center of the screen.  This amounts to figuring out how much to pan in the up and right directions. 

If you play with the camera.pan() function you will notice that it appears to use screen coordinates as its arguments.  This makes sense from a GUI interface standpoint in that when the user pans the view, they typically do so with the mouse and therefore start and end screen coordinates of a mouse movement are the most natural inputs to a pan routine.  Unfortunately for us, that makes panning more of a challenge for us.  Nonetheless, lets approach the problem from an iterative viewpoint similar to what we did with the rotations. 

We’re going to attach the panning problem using a similar approach as we used for the rotate view function.  That is, we’re going to start with a trial pan in one of the directions, look to see what that does in world coordinates, then perform a real pan in that direction, or its direct opposite, to try to get closer to our target.  In a similar fashion to the rotate code, we project the pan direction vector in world coordinates onto a plane defined by our current view vector.  At each iteration we test the vector from our focus point to the desired point to see if it is parallel to the view vector.  If it is, then that means our desired point lies right down the axis of the view vector passing through the focus point.  Thus, it would be in the center of the screen.  The following code implements this technique.

function pan_view(desired_point) {
    // Get the camera and graphics objects
    var camera = DS.Graphics.Camera;
    var graphics = DS.Graphics;
    // Maximum amount to pan by
    var height = DS.Graphics.GfxWindow.height;
    var width = DS.Graphics.GfxWindow.width;
    // The current focus point
    var fp = null;
    // The trial focus point
    var fp_trial = null;
    // The trial amount to pan
    var trial_pan = (height + width) / 40;  //Unit is pixels...
    // View direction
    var view_v = null;
    // Vector from the focus point to the desired point
    var fp_to_dp = null;
    // Vector of pan direction
    var pan_dir = null;
    // The amount we are going to pan
    var pan_amount = null;
    // Tolerance for determining if the view vector is parallel
    // to the focus point to desired point vector
    var eps = 1e-4;
    // We'll flip back and forth between panning up and panning right
    var b_first_arg = true;
    // The current count of our iteration
    var cnt = 0;
    // The maximum number of iterations to perform
    var max_iter = 200;

    var look_at = desired_point;
    var max_pan = (height + width) / 20;
    do {
        fp = [camera.FocusPointX,
              camera.FocusPointY,
              camera.FocusPointZ];
        vd = normalize([camera.ViewDirectionX,
                        camera.ViewDirectionY,
                        camera.ViewDirectionZ]);

        if (b_first_arg == true) {
            camera.pan(trial_pan, 0);
        } else {
            camera.pan(0, trial_pan);
        }
        trial = [camera.FocusPointX,
                 camera.FocusPointY,
                 camera.FocusPointZ];
        // In world coordinates, see which direction this pan took us
        pd = sub(trial, fp);
        // Project this onto a plane defined by the current view 
        // vector
        pd_in_view = cross(vd, cross(pd, vd));

        // Move back to our original location
        if (b_first_arg == true) {
            camera.pan(-trial_pan, 0);
        } else {
            camera.pan(0, -trial_pan);
        }
        fp = [camera.FocusPointX,
              camera.FocusPointY,
              camera.FocusPointZ];
        // Take a vector from the current focus point to the
        // desired center point
        fp_to_dp = normalize(sub(look_at, fp));
        // Only try to pan the view if the vector between our focus
        // point and view vector are not parallel.
        if (Math.abs(1.0 - Math.abs(dot(fp_to_dp, vd))) > eps) {
            // Project this onto a plane defined by the view direction
            fp_to_dp_in_view = cross(vd, cross(fp_to_dp, vd));

            // Normalize the pan test vector so we can walk our way there
            pd_in_view = normalize(pd_in_view);

            // Figure out how much we need to pan

            pan_amount = trial_pan * dot(fp_to_dp_in_view, pd_in_view);
            // Knock it down by a factor associated with our iteration
            // scheme.  This helps prevent spurious jittering as we 
            // make our way there.
            var factor = cnt / max_iter;
            pan_amount *= (1.0 - factor * factor * factor);

            if (b_first_arg == true) {
                camera.pan(pan_amount, 0);
            } else {
                camera.pan(0, pan_amount);
            }
            b_pan_converged = false;
        } else {
            b_pan_converged = true;
        }
        // Flip the direction
        b_first_arg = !b_first_arg;
//        graphics.Refresh();
    } while (cnt++ < max_iter && !b_pan_converged);
}

If you’ve hung with me this far, you know that we’ve got a function that will rotate our view so that we’re looking down a particular axis.  We’ve also got a function that will rotate the view about our line of sight so that we can make sure a particular direction is up.  Finally, we’ve got a function that will pan the view around so that a particular point in space shows up right in the middle of the screen.  All of this was accomplished with the limited set of functions and properties we could deduce from the Camera object using the debugger.  We also demonstrated how to create a slice plane using the technique of “drawing” al line straight across the middle of the screen.  With these four functions we should have all we need conceptually to create a slice plane right where we want it, oriented exactly as we want it. 

Implementation

The following code implements these ideas and creates a little dialog box that is presented to the user.  I’ll cover the implementation of a dialog box in a separate post, however, you may be able to decode the secret sauce here anyway.   The user interface works as follows.  You pick a coordinate system from the drop down list for which you would like to create a slice plane.  You can then either create a slice plane from this coordinate system, or have the code automatically reorient the view so that you are looking down the Z axis of this coordinate system.   Below is a screen shot of the dialog box that is presented, and finally, a listing of the entire code for the script.

image

Listing of the Code

Here is a listing of the entire script.

// This string contains the HTML code for the dialog box
// we will present to the user
var html_text = 
'<html xmlns="http://www.w3.org/1999/xhtml" >\n'+
'<head>\n'+
'<title>Slice Plane Creation Tool</title>\n'+
'<style type="text/css">\n'+
'    body\n'+
'    {\n'+
'        background: buttonface;\n'+
'        font: messagebox;\n'+
'    }\n'+
'    table\n'+
'    {\n'+
'        border-collapse:collapse;\n'+
'        font: messagebox;\n'+
'    }\n'+
'    select\n'+
'    {\n'+
'        margin: 2px 8px 2px 4px;\n'+
'        width: 250px;    \n'+
'        height: 28px;\n'+
'    }\n'+
'    input\n'+
'    {\n'+
'        margin: 0px 4px 0px 0px;\n'+
'    }\n'+
'</style>\n'+
'<script type="text/javascript">\n'+
'    var wb, xfer;\n'+
'    function window_onload() {\n'+
'        xfer = window.external;\n'+
'        wb = xfer("wb");\n'+
'        var cs_group = xfer(\'csys_group\');\n'+
'        var cs_select = xfer(\'csys_select\');\n'+
'        if (cs_group != null) {\n'+
'            populate_coordinate_systems(cs_group, cs_select);\n'+
'        }\n'+
'    }\n'+
'    function window_onunload() {\n'+
'        xfer("CSYS") = parseInt(csys_select.value);\n'+
'        xfer("CreateSlicePlane") = create_plane.checked;\n'+
'        xfer("ViewPlane") = view_plane.checked;\n'+
'    }\n'+
'    function populate_coordinate_systems(cs_group, cs_select) {\n'+
'        csys_select.options.length=0;\n'+
'        for (var i = 1; i <= cs_group.Children.Count; i++) {\n'+
'            var csys = cs_group.Children.Item(i);\n'+
'            var csys_name = csys.Name;\n'+
'            var csys_id = csys.ID;\n'+
'            if (csys_id == cs_select) {\n'+
'                csys_select.options[i - 1] = new Option(csys_name, csys_id, true, true);\n'+
'            } else {\n'+
'                csys_select.options[i - 1] = new Option(csys_name, csys_id, false, false);\n'+
'            }\n'+
'        }\n'+
'    }\n'+
'</script>\n'+
'</head>\n'+
'<body onload="window_onload()" onunload="window_onunload()" scroll="no">\n'+
'<table>\n'+
'<tr>\n'+
'<td>\n'+
'<table>\n'+
'<tr>\n'+
'<td nowrap="true">Coordinate System: </td>\n'+
'<td><select id="csys_select"></select></td>\n'+
'</tr>\n'+
'</table>\n'+
'</td>\n'+
'</tr>\n'+
'<tr>\n'+
'<td>\n'+
'<input id="create_plane" type="checkbox" checked="checked" />Create Slice Plane\n'+
'</td>\n'+
'</tr>\n'+
'<tr>\n'+
'<td>\n'+
'<input id="view_plane" type="checkbox"  />Look at Plane\n'+
'</td>\n'+
'</tr>\n'+
'</table>\n'+
'</body>\n'+
'</html>\n';


var SC = DS.Script;

// Entry point for the whole script.
function main() {
 
    // Create the dialog
    var slice_plane_dialog =
        SC.CreateActiveXObject(
            SC.GenWBProgId('WBControls.WBHTMLDialog')
            );
    var dlg_OK = 1;
    var dlg_Cancel = 2;
    var dlg_VScroll = 4;
    var flags = dlg_OK | dlg_Cancel | dlg_VScroll;
    var caption = 'Slice Plane From Coordinate System';
    var width = 410;
    var height = 110;
    var b_modal = false;
    // Create the html code from a temporary file
    var fso = SC.fso;
    var tfolder = fso.GetSpecialFolder(2);
    var tname = fso.GetTempName();
    var tfile = tfolder.CreateTextFile(tname, true);
    tfile.WriteLine(html_text);
    tfile.Close();
    var path = fso.BuildPath(tfolder.Path, tname);
    var xfer =
        SC.CreateActiveXObject(
            SC.GenWBProgId('WBControls.DlgArgs')
            );

    xfer('wb') = WB;
    xfer('dlg') = slice_plane_dialog;
    var csys_group = get_coordinate_system_group();
    
    if(csys_group == null) {
        SC.WBScript.Out('Cannot find the coordinate system group in the tree.' +
                        ' Please create a coordinate system first', true);
        return;
    }
    // Pass over to the dialog the coordinate system group.
    xfer('csys_group') = csys_group;
    var active_obj = DS.Tree.FirstActiveObject;
    if (active_obj.Class == SC.id_CoordinateSystem) {
        xfer('csys_select') = active_obj.ID;
    } else {
        xfer('csys_select') = -1;
    }
    // Show the dialog
    var ret = slice_plane_dialog.DoDialog(WB.hWnd,
                                          path,
                                          xfer,
                                          b_modal,
                                          width,
                                          height,
                                          flags, 
                                          caption);
    if (ret == 1) {
        // We get here if the user presses OK on the dialog, 
        // so lets grab the user's choices back from the selection.
        var selected_csys = xfer('CSYS');
        var create_plane = xfer('CreateSlicePlane');
        var look_at= xfer('ViewPlane');
        // Do the work here
        var origin = get_coordinate_system_origin(selected_csys);
        var z_axis = get_coordinate_system_z_axis(selected_csys);
        var y_axis = get_coordinate_system_y_axis(selected_csys);
        // Bail out if we can't create the vectors for the coordinate system
        if (origin == null || z_axis == null || y_axis == null) {
            SC.WBScript.Out('Cannot determine the coordinate system vectors.' + 
                            ' Unable to create slice plane', true);
            return;
        }
        if (create_plane == true) {
            create_slice_plane(origin, z_axis, y_axis);
        }
        if (look_at == true) {
            look_at_plane(origin, z_axis, y_axis);
        }
    }

    // Delete the html file to clean up
    fso.DeleteFile(path);
}

// Very simple utility function that is used to write out an html
// file to another file in javascript string format.  This is useful
// for pasting into a source file to package the html file into
// a single source script for distribution purposes.
function write_html(path) {
    var fso = SC.fso;
    var input = fso.OpenTextFile(path, 1, false);
    var output = fso.CreateTextFile(path + '.out', true);

    while (!input.AtEndOfStream) {
        var line = input.ReadLine();
        output.WriteLine('\'' + line + '\\' + 'n\'+');

    }
    output.WriteLine('\'\';');
    output.Close();
    input.Close();
}    

// Return the origin of a given coordinate system
function get_coordinate_system_origin(csys_id) {
    var csys = get_coordinate_system(csys_id);
    if (csys == null) {
        return null;
    }
    return [csys.OriginXLocation, 
            csys.OriginYLocation,
            csys.OriginZLocation];
}

// Return a vector oriented along the z axis of a given
// coordinate system
function get_coordinate_system_z_axis(csys_id) {
    var csys = get_coordinate_system(csys_id);
    if (csys == null) {
        return null;
    }
    return normalize([csys.ZDirectionXValue,
                      csys.ZDirectionYValue,
                      csys.ZDirectionZValue]);
}

// Return a vector oriented along the y axis of a given
// coordinate system
function get_coordinate_system_y_axis(csys_id) {
    var csys = get_coordinate_system(csys_id);
    if (csys == null) {
        return null;
    }
    return normalize([csys.YDirectionXValue,
                      csys.YDirectionYValue,
                      csys.YDirectionZValue]);
}

// Return an actual coordinate system object given
// a tree ID.
function get_coordinate_system(csys_id) {
    var group = get_coordinate_system_group();
    if (group == null) {
        return null;
    }
    for (var i = 1; i <= group.Children.Count; i++) {
        var child = group.Children.Item(i);
        if (child.ID == csys_id) {
            return child;
        }
    }
    return null;
}

// Return the coordinate system group in the tree.
function get_coordinate_system_group() {
    var branch = SC.getActiveBranch();
    if (branch == null) {
        return null;
    }
    return branch.CoordinateSystemGroup;
}

// Create a slice plane given a point, plane normal and up
// vector.  This function uses the orientation functions
// below
function create_slice_plane(point, plane_normal, inplane_up) {
    // Get a handle to the slicetool
    var slice_tool = DS.Graphics.SliceTool;

    // Orient the view using our view orientation functions
    for (var i = 0; i < 3; i++) {
        rotate_view(inplane_up);
        pan_view(point);
        rotate_up(plane_normal);
        DS.Graphics.Refresh();
    }

    // Get the window's dimensions
    var height = DS.Graphics.GfxWindow.height;
    var width = DS.Graphics.GfxWindow.width;

    // Create the slice plane
    var slice_plane =
        slice_tool.SetPlane(0, height / 2, width, height / 2);
    DS.Graphics.Refresh();
}

// Look at a particular point along an axis.  This is useful
// for looking at a slice plane after it has been created.
function look_at_plane(point, normal, up) {
    // Orient the view using our view orientation functions

    // Orient the view using our view orientation functions
    for (var i = 0; i < 3; i++) {
        rotate_view(normal);
        pan_view(point);
        rotate_up(up);
        DS.Graphics.Refresh();
    }

}

// Calculate the dot product of two vectors
function dot(va, vb) {
    return va[0] * vb[0] + va[1] * vb[1] + va[2] * vb[2];
}

// Calculate the length of a vector
function norm(va) {
    return Math.sqrt(va[0] * va[0] + va[1] * va[1] + va[2] * va[2]);
}

// Normalize a vector
function normalize(va) {
    var len = norm(va);
    if (len != 0.0) {
        return [va[0] / len, va[1] / len, va[2] / len];
    }
}

// Calculate the cross product of two vectors
function cross(va, vb) {
    return [va[1] * vb[2] - va[2] * vb[1],
            va[2] * vb[0] - va[0] * vb[2], 
            va[0] * vb[1] - va[1] * vb[0]];
}

// Scale a vector by a scalar
function scale(va, sc) {
    return [va[0] * sc, va[1] * sc, va[2] * sc];
}

// Add two vectors
function add(va, vb) {
    return [va[0] + vb[0],
            va[1] + vb[1],
            va[2] + vb[2]];
}

// Subtract vector b from vector a
function sub(va, vb) {
    return [va[0] - vb[0],
            va[1] - vb[1],
            va[2] - vb[2]];
}

function are_parallel(va, vb, tol) {
    return 1.0 - (Math.abs(dot(va, vb)) / (norm(va) * norm(vb))) < tol ? true : false;
}

function are_perpendicular(va, vb, tol) {
    return (Math.abs(dot(va, vb)) / (norm(va) * norm(vb))) < tol ? true : false;
}

// Rotate the view around so that we are looking in the 
// direction of the desired view
function rotate_view(desired_view) {
    // Get the camera and graphics objects
    var camera = DS.Graphics.Camera;
    var graphics = DS.Graphics;
    // This is our tolerance for being parallel to the desired view
    var eps = 1e-7;
    // This is the maximum number of iterations we'll try.  
    // While loops with a tolerance check only scare me...
    var max_iter = 200;
    // Get the current view
    var view_c = normalize([camera.ViewDirectionX,
                            camera.ViewDirectionY,
                            camera.ViewDirectionZ]);
    // Make sure we're normalized
    dvd = normalize(desired_view);
    // This should be close to 1 if parallel.

    var cnt = 1;
    var b_first_arg = true;
    var normal = null;
    var trial = null;
    var view_p = null;
    var trial_ang = 15;
    var applied_rotation = null;
    var previous_up = 0;
    var previous_right = 0;
    var right_factor = up_factor = -1.0;
    // Loop until we're parallel, or we give up trying
    do {
        var factor = cnt / max_iter;
        factor *= factor;
        factor *= factor;
        // Get the view direction
        vd = normalize([camera.ViewDirectionX,
                        camera.ViewDirectionY,
                        camera.ViewDirectionZ]);
        // Get the current up vector
        camera.rotate(trial_ang, 0);
        trial = normalize([camera.ViewDirectionX,
                       camera.ViewDirectionY,
                       camera.ViewDirectionZ]);
        current_up = normalize(cross(trial, vd));
        // Rotate back so we don't lose our place
        camera.rotate(-trial_ang, 0);
        

        vd = normalize([camera.ViewDirectionX,
                        camera.ViewDirectionY,
                        camera.ViewDirectionZ]);

        if (b_first_arg) {
            camera.rotate(trial_ang, 0);
        } else {
            camera.rotate(0, trial_ang);
        }
        trial = normalize([camera.ViewDirectionX,
                           camera.ViewDirectionY,
                           camera.ViewDirectionZ]);
        n = normalize(cross(vd, trial));
        // Rotate back
        if (b_first_arg) {
            camera.rotate(-trial_ang, 0);
        } else {
            camera.rotate(0, -trial_ang);
        }
        vd = normalize([camera.ViewDirectionX,
                        camera.ViewDirectionY,
                        camera.ViewDirectionZ]);
        // Only do this rotation if our desired view vector
        // is not nearly parallel to our current axis of rotation
        // If it is, then the code inside the if statement will
        // be skipped.
        if (!are_parallel(dvd, n, eps)) {
            // Project the desired view vector onto our plane of rotation
            vp = normalize(cross(n, normalize(cross(dvd, n))));
            // Calculate the angle between the projected vector
            // and our current view direction
            dot_product = dot(vp, vd);
            needed_rotation = Math.acos(dot(vp, vd)) * 180 / Math.PI
            
            // Knock it down by a factor associated with our iteration
            // scheme.  This helps prevent spurious jittering as we 
            // make our way there.
            needed_rotation *= (1.0 - factor);
            if (b_first_arg) {
                // If we start to diverge, try rotating back the other way
                if(previous_up < needed_rotation) {
                    up_factor = (up_factor < 0.0) ? 1.0 : -1.0;
                    previous_up = needed_rotation;
                }
                camera.rotate(up_factor * needed_rotation, 0.0);
                
                
            } else {
                if (previous_right < needed_rotation) {
                    right_factor = (right_factor < 0.0) ? 1.0 : -1.0;
                    previous_right = needed_rotation;
                }
                camera.rotate(0.0, right_factor * needed_rotation);
            }
        }
        // See if we are there yet
        vd = normalize([camera.ViewDirectionX,
                        camera.ViewDirectionY,
                        camera.ViewDirectionZ]);
        b_rotate_converged = are_parallel(vd, dvd, eps);
        if (b_rotate_converged) break;
        // Flip directions
        b_first_arg = !b_first_arg;
//        DS.Graphics.Refresh();
    } while (cnt++ < max_iter);
        
}

function rotate_up(desired_up_vector) {
    
    // Get the camera and graphics objects
    var camera = DS.Graphics.Camera;
    var graphics = DS.Graphics;
    var view_c = null;
    var trial = null;
    var current_up = null;
    var theta = null;
    var trial_ang = 15;
    var theta1;
    var theta2;
    var eps = 0.01;
    // Make sure our passed in value is normalized
    duv = normalize(desired_up_vector);

    // Which way are we currently looking now?
    vd = normalize([camera.ViewDirectionX,
                    camera.ViewDirectionY,
                    camera.ViewDirectionZ]);

    // Perform our trial rotation about the screen vertical axis
    camera.rotate(trial_ang, 0);
    trial = normalize([camera.ViewDirectionX,
                       camera.ViewDirectionY,
                       camera.ViewDirectionZ]);
    up = normalize(cross(trial, vd));
    // Rotate back so we don't lose our place
    camera.rotate(-trial_ang, 0);
    // How much to we need to rotate?
    theta1 = Math.acos(dot(duv, up)) * 180 / Math.PI;
    // Rotate assuming a positive rotation
    camera.rotateZ(theta1);
    DS.Graphics.Refresh();
    // Now we need to test to see if this is indeed the proper direction
    // Which way are we currently looking now?
    vd = normalize([camera.ViewDirectionX,
                    camera.ViewDirectionY,
                    camera.ViewDirectionZ]);
    // Perform our trial rotation about the screen vertical axis
    camera.rotate(trial_ang, 0);
    trial = normalize([camera.ViewDirectionX,
                       camera.ViewDirectionY,
                       camera.ViewDirectionZ]);
    up = normalize(cross(trial, vd));
    // Rotate back so we don't lose our place
    camera.rotate(-trial_ang, 0);
    // How much to we need to rotate?
    theta2 = Math.acos(dot(duv, up)) * 180 / Math.PI;
    if (theta2 > theta1) {
        // Rotate backwards if our second guess is larger than our first
        camera.rotateZ(-theta2);
//        DS.Graphics.Refresh();
    }
}

function pan_view(desired_point) {
    // Get the camera and graphics objects
    var camera = DS.Graphics.Camera;
    var graphics = DS.Graphics;
    // Maximum amount to pan by
    var height = DS.Graphics.GfxWindow.height;
    var width = DS.Graphics.GfxWindow.width;
    // The current focus point
    var fp = null;
    // The trial focus point
    var fp_trial = null;
    // The trial amount to pan
    var trial_pan = (height + width) / 40;  //Unit is pixels...
    // View direction
    var view_v = null;
    // Vector from the focus point to the desired point
    var fp_to_dp = null;
    // Vector of pan direction
    var pan_dir = null;
    // The amount we are going to pan
    var pan_amount = null;
    // Tolerance for determining if the view vector is parallel
    // to the focus point to desired point vector
    var eps = 1e-4;
    // We'll flip back and forth between panning up and panning right
    var b_first_arg = true;
    // The current count of our iteration
    var cnt = 0;
    // The maximum number of iterations to perform
    var max_iter = 200;

    var look_at = desired_point;
    var max_pan = (height + width) / 20;
    do {
        fp = [camera.FocusPointX,
              camera.FocusPointY,
              camera.FocusPointZ];
        vd = normalize([camera.ViewDirectionX,
                        camera.ViewDirectionY,
                        camera.ViewDirectionZ]);

        if (b_first_arg == true) {
            camera.pan(trial_pan, 0);
        } else {
            camera.pan(0, trial_pan);
        }
        trial = [camera.FocusPointX,
                 camera.FocusPointY,
                 camera.FocusPointZ];
        // In world coordinates, see which direction this pan took us
        pd = sub(trial, fp);
        // Project this onto a plane defined by the current view 
        // vector
        pd_in_view = cross(vd, cross(pd, vd));

        // Move back to our original location
        if (b_first_arg == true) {
            camera.pan(-trial_pan, 0);
        } else {
            camera.pan(0, -trial_pan);
        }
        fp = [camera.FocusPointX,
              camera.FocusPointY,
              camera.FocusPointZ];
        // Take a vector from the current focus point to the
        // desired center point
        fp_to_dp = normalize(sub(look_at, fp));
        // Only try to pan the view if the vector between our focus
        // point and view vector are not parallel.
        if (Math.abs(1.0 - Math.abs(dot(fp_to_dp, vd))) > eps) {
            // Project this onto a plane defined by the view direction
            fp_to_dp_in_view = cross(vd, cross(fp_to_dp, vd));

            // Normalize the pan test vector so we can walk our way there
            pd_in_view = normalize(pd_in_view);

            // Figure out how much we need to pan

            pan_amount = trial_pan * dot(fp_to_dp_in_view, pd_in_view);
            // Knock it down by a factor associated with our iteration
            // scheme.  This helps prevent spurious jittering as we 
            // make our way there.
            var factor = cnt / max_iter;
            pan_amount *= (1.0 - factor * factor * factor);

            if (b_first_arg == true) {
                camera.pan(pan_amount, 0);
            } else {
                camera.pan(0, pan_amount);
            }
            b_pan_converged = false;
        } else {
            b_pan_converged = true;
        }
        // Flip the direction
        b_first_arg = !b_first_arg;
//        graphics.Refresh();
    } while (cnt++ < max_iter && !b_pan_converged);
}

main();

ANSYS Mechanical Scripting: HOWTO Part 4

Where we left Off Last Time…

In the third part of this series, I left you with an introduction to the ANSYS Selection Manager and Geometry B-Rep.  We also created a sample script that iterated over all of the faces in a selected body and created a named selection for each face.   We built on our technique described in Part 2 for discovering ANSYS functionality. 

Introducing the Microsoft Visual Studio Script Debugger

In this post I want to introduce you to the Visual Studio Script Debugger.  Running ANSYS Mechanical underneath the debugger is definitely the most powerful tool we have at our disposal when it comes to creating custom scripts.  If you can learn to use this tool effectively, you will be able to watch the inner workings of Mechanical interact with your scripts in a precisely controlled environment.  Furthermore, the rich debugger displays will allow you to see most all of the mechanical objects in suspended animation and will also allow you to investigate their properties and methods in real time.  Finally, the “Immediate Window” in the debugger will allow you to execute random pieces of code that you type in on the fly while Mechanical is suspended in the debugger.  This will allow you to test out various techniques and code sequences in one debugging session.   Simply put, you must become proficient at using the debugger to have any real hope of successfully scripting the ANSYS Mechanical tool.

Some Preliminaries

There are a number of steps you need to take to actually get the debugger environment set up properly.  Here are a list of steps you need to take:

  1. You obviously need to install Visual Studio.  I use Visual Studio Professional 2010, but I think any version including 2005 and above will probably be sufficient.  I cannot speak to whether or not the Express Editions will work for this task, so if someone has experience with one of these editions and can shed some light, please post a comment.
  2. In the Internet Explorer “Internet Options” dialog (arrived at through Tools->Internet Options…) you need to make sure that “Disable Script Debugging (Other)” is unchecked.  See image below.
  3. Finally, you need to adjust the registry setting HKEY_CURRENT_USER\Software\Microsoft\Windows Script\Settings\JITDebug = 1

 

image

 

Debugging Code

The easiest way to stop Mechanical in the debugger is to insert the javascript “debugger” command wherever you would like to stop execution and look around.  So, for example, lets assume we have the following script:

function my_func() {
    debugger;
    var my_array = new Array();
    for(var i = 0; i < 10; i++) {
        my_array[i] = 2*i;
    }
}

my_func();

As you can see above, this simple script consists of one function called my_func that just populates an array.  The first statement in this function is the keyword debugger.   If all works according to plan, when we execute this script, we’ll end up in the Microsoft Debugger.  To make this all happen, you need to follow a ten step program:

  1. Start MS Visual Studio
  2. Start ANSYS Mechanical
  3. Copy the above script to a file named debug_test.js anywhere on your harddrive.
  4. Inside MS Visual Studio, navigate to Debug->Attach to Process…image
  5. In the dialog that pops up, look for the AnsysWBU.exe executable in the list.  image
  6. Select AnsysWBU.exe and press the Attach button.
  7. Switch to the Mechanical Application
  8. From the tools Menu, choose Tools->Run Macro…
  9. Navigate to the location of debug_test.js and click open.
  10. Start debugging…

At this point you can set breakpoints, step into and out of functions, interrogate variable values, etc…  The debugger is truly a remarkable tool, as we shall briefly touch on next.

Whirlwind Tour of the MS Debugger

Below is a screen shot of our little script inside the visual studio debugger.

image

As you can see, the current execution line is located over the debugger statement.  Pressing F10 will advance one line at a time.  Also, you will notice that there is a window below the code window with a tab called “Locals”  This tab shows the local variables defined in this function.  At this point in the execution, they are not defined.  However, as we step the execution, these variables will come into play.  You can watch what happens to variables as you step through your code here.

Also, you will notice a second window called the “Immediate Window”.  This window may not be visible for your default debugging layout.  If you don’t see it, navigate to Debug->Windows->Immediate Window

image

In the immediate window you can enter in code that will be executed “immediately” while the program is stopped in the debugger.  So, for example, if you see a problem in your script that might be fixed with a code change you may be tempted to leave the debugger fix it and then re-run.  However, if you are unsure of your fix, or if you want to try a number of alternatives, you can just type them into the immediate window and see what happens.  It works as if you were executing new code at the current breakpoint location that isn’t in your script.  You can then step forward in your script to execute a few lines, enter something in the immediate window, and then continue on stepping.  It is a very powerful technique.

Lets step through a few lines.  I’m going to press F10 four or five times.  At that point I get the following screen.

image

You can see that my_array now has an ellipse next to it and I is equal to 1.  If I click on the plus symbol next to my_array in the locals pane, I get the following:

image

You can see that we’ve filled the zero-th item in the array with the value zero.  Lets hit F10 some more…  Note, you can also hover your mouse over the array in the code window.  You can see I’ve filled in a few more values

image

The print screen function doesn’t capture it well, but you can see the values in the array have changed.  Now, lets say I just realized that for some code below to work properly, the fifth value in this array needs to be less than 10.  (We don’t have any code below, and this is quite contrived, but bear with me and use your imagination…)  So, I just pop down into the immediate window and type:

my_array[5] = 8;

image

When I press return, this will execute that single line of code. You can see that my_array[5] now holds the value 8 instead of 10.

image

 

You can also set break points in your code.    Just right click on the line, and choose Breakpoint->Insert Breakpoint

image

Breakpoints work kind of like the debugger command, only they don’t litter your code with commands.  I usually insert one debugger command at the top of my script and then breakpoints where needed.   Breakpoints can also be conditional.  That is, they won’t break unless some condition is true.  This is great for executing most of a loop, but breaking on the last or next to last item.  After you have inserted a breakpoint, just right click on the line and choose Breakpoint->Condition…

image

You can insert a condition like:

image

From here, press F5 which is the Run command, and the code will execute inside the loop until I == 8.

Conclusion

The debugger is your friend; plain and simple.  Use it wisely and use it often.  In fact, if you want to script in Mechanical, you are going to have to become intimately familiar with the debugger.  In a future article, I will show you more powerful debugging tips when the debugger by itself lets you down. Stay tuned…

ANSYS Mechanical Scripting: HOWTO Part 3

Where We Left Off Last Time…

In the second part of this series I left you with a recipe for discovering ANSYS Mechanical functionality by using entities in the GUI itself to help guide your search through the ANSYS source code that is shipped with every ANSYS Mechanical installation.  This technique is the first tool in our toolbox.  In this blog post, I’ll introduce you to a couple of the internal ANSYS Mechanical data structures that are associated with geometry in particular.  I’ll simply be introducing them in this blog without any discussion of how one might determine additional functionality associated with these structures.  In a future blog post, I’ll step you through the use of the most powerful tool available in our toolbox: the Visual Studio Script Debugger.

Introducing the ANSYS Mechanical Geometry Selection Manager and B-Rep Geometry Data Structure

CAD geometry is represented inside a computer using a Boundary Representation (B-Rep) data structure.  This data structure falls into a category of data structures computer scientists like to call graphs.   In general, a graph is a structure comprised of nodes connected to each other with arcs.  A road map is a type of graph where cities represent nodes and the roads connecting them represent arcs.  Sites like facebook and linkedin are giant graphs where people are represented as nodes and friendships or acquaintances are represented as arcs.   As you can see, graphs are pretty general and you can represent lots of structured relationships with them.

In CAD geometry, the nodes in our graph are comprised of various topological entities like vertices, edges, faces, shells and bodies.  The arcs are the incident relationships that exist between these entities.  So, for example, an edge will have a start vertex and an end vertex.  These relationships can be represented by links between the particular edge and its start vertex or end vertex.  Likewise, however, the particular start vertex of one edge may be the end vertex of another edge.  So, this relationship can be represented by links between a given vertex and all of the edges to which it is attached.  There are similar relationships between edges and faces.  A series of edges will form the boundary for a given face and a series of faces might share a common edge.  Shells represent a volume that is defined by a set of connected faces. (i.e faces that share common edges)  A body is comprised of one or more shells.  Why isn’t a body the same as a shell?  Well, it turns out that bodies with internal holes are more conveniently represented as comprising of two or more shells.  One shell represents the outermost envelope and the other shells represent the internal voids.   How, then, does one know if they are inside or outside of a body?  Well, the faces that comprise a shell each have an associated outward normal vector that points to the “outside” of the body in question.  The normal vector sense is defined using an ordering on the bounding edges of a given face.  Typically as you traverse the list of edges comprising the boundary of a given face, you are actually walking counter-clockwise around the face.  Using this assumption, one can determine an outward normal.

Lets discuss some of the specifics of the ANSYS Mechanical B-Rep.  First, an individual B-Rep structure is associated with each part in a Mechanical session.  So, if you have 15 parts you will have 15 different B-Rep structures.  Second, each entity (vertex, edge, face, body) within a given B-Rep has an associated unique ID. (This ID is called the TopoID, or Topological ID)  This ID is just an unsigned integer value, with one small twist.  The top two bits of the ID encode the entity’s topological type.  That is, if you stripped off the top two bits, you would get one of the following four combinations:

Top Two Bits

Decimal Value

Topological Type

00 0 Vertex
01 1 Edge
10 2 Face
11 3 Body

The lower 30 bits actually encode the unique id for a given vertex, edge, face or body.  So, what that means is that ANSYS Mechanical is limited to geometry with fewer than roughly one billion vertices, edges, etc… If you have a model with more than that, perhaps a bit of geometry simplification is in order.

To access a given topological entity within a B-Rep object, one uses the following type of function call (Note, I haven’t introduced where the brep_object below comes from…  Bear with me, we’ll get to it soon.)

var topo_object = brep_object.Cell(topo_id);

The Cell function accepts a valid topoID and returns the associated object.  The returned object is a rich object and has numerous fields that obviously differ depending on the topological type.  So, for example, an edge object will have a StartVertex and EndVertex field that contains references to its start and end vertices in particular.  However, all topological objects have one field in common.  They all contain an ID field.  So the following code would supply the topolD of the start vertex associated with a given edge:

var brep_object = /*Get B-Rep from somewhere */
var edge_id = /* Get Edge ID from somewhere*/
var edge = brep_object.Cell(edge_id);
var start_vertex_id = edge.StartVertex.Id;

You might gather that one can get the actual start vertex object with a subsequent line like:

var start_vertex = brep_object.Cell(start_vertex_id);

Now, each vertex object is going to contain a list of edges that are attached to this vertex.  So, if we needed all the edges attached to the original edge for some reason, we could use this technique with the start vertex and the end vertex. 

Walking the B-Rep data structure is perhaps interesting at best, but it currently doesn’t do anything for us.  What we need is to be able to select some piece of geometry, or pieces of geometry.  Enter the ANSYS Mechanical Selection Manager.  This object is closely associated with the B-Rep data structures.  Before we discuss this object in some detail, consider starting each and every one of your scripts with the following code:

//  The script code object is a reference to the javascript engine embedded
//  inside the mechanical application.  This object has all of the functions
//  that are defined in the ANSYS mechanical javascript files
var SC = DS.Script;

//  The selection manager object is a reference to the ANSYS mechanical 
//  geometry selection tool.  This is a complex object that contains 
//  references to the currently selected set of geometric entities.  More
//  importantly, it contains back pointers to entity ids within the 
//  Mechanical BREP structure.
var SM = SC.sm;

This code introduces two global variables into your script.  The first object contains a reference to the actual javascript engine embedded in the Mechanical application.  This object is derived from the DS object, which represents the fundamental Mechanical application object.  You can think of the DS object as the mother of all objects.  The ScriptCode or SC object stores all of the javascript code that ships with the ANSYS mechanical product.  So, if you use the technique described in Part 2 of this series and find a particular function call, you would actually call it against this object.  So, for example, we found the function doPrototypeInsertCommandEditor() in Part 2.  If we were to actually call that function in our script, it would look like:

SC.doPrototypeInsertCommandEditor();

The second global variable shown above is the selection manager, or SM.  You can see that this object is a child of the script code object.   This object has numerous functions and properties and in a later post I hope to document as many as I can.  However, a great deal of mileage can be gained by understanding the following two function calls, and the following few properties.  I list them below in code form so you can see how they would typically be called inside a script:

/*  F U N C T I O N S   O F   T H E    S E L E C T I O N    M A N A G E R */
// Clear()  
// This function resets the selection manager so that no geometry is currently 
// selected.
SM.Clear();

// AddToSelection(partID, topoID, bShouldDraw)
// This function adds a particular topological entity (vertex, edge, face, body) 
// to the currently selected set. The first argument specifies the particular
// part this entity is associated with.  You can think of this as referencing the
// particular BRep structure.  The second argument specifies the topological
// entities specific ID.  The final argument is a boolean specifying whether or
// not the selection manager should highlight the selected geometry in green.
SM.AddToSelection(partID, topoID, bShouldDraw);

/*  P R O P E R T I E S   O F   T H E    S E L E C T I O N    M A N A G E R */
// Vertices
// The number of currently selected vertices
var num_selected_vertices = SM.Vertices;

// Edges
// The number of currently selected edges.
var num_selected_edges = SM.Edges;

// Faces
// The number of currently selected faces.
var num_selected_faces = SM.Faces

// Parts
// The number of currently selected (bodies*)
// Note: I think this is where the SM and the BRep use different
// terms for a body.  However, I may be wrong.
var num_selected_bodies = SM.Parts

// SelectedCount
// The total number of selected entities.  (Sum of the above)
var num_selected = SM.SelectedCount;

// SelectedEntityTopoID(i)
// The ith selected entity's topological id.  Note you can use the
// top two bits of this id to determine the type of this selected
// object.  i goes from 1 to SelectedCount.
var selected_topoid = SM.SelectedEntityTopoID(i);

// SelectedPartID(i)
// The part id for the ith selected entity. For models with numerous
// parts this will tell you which part the currently selected entity is
// associated with.  You use this to reference the proper BRep structure.
// i goes from 1 to SelectedCount.
var selected_partid = SM.SelectedPartID(i);

// PartMgr
// The Part Manager (PartMgr) object is a child of the selection manager.
// It also has numerous methods and properties, however, the one
// property that I have found most useful is the PartById(id).
// PartMgr.PartById(id)
// This returns the particular part object associated with a given part id
var part = SM.PartMgr.PartById(partID);

// The part object returned by the PartMgr.PartById(partID) call has 
// associated with it a given BRep data structure.  Referencing this brep
// structure is how you get from the selection manager's selected entity
// id to a given piece of topology.  So for example, assume the user
// has selected one face in the model.  This is one way we can get to that
// face object.
var topo_id = SM.SelectedEntityTopoID(1);
var part_id = SM.SelectedPartID(1);
var part = SM.PartMgr.PartById(part_id);
var brep = part.BRep;
var face = brep.Cell(topo_id);

The above code represents a small portion of the interface of the ANSYS Mechanical selection manager and also the geometry B-Rep structure.  As you can see, the two structures are coupled through the notion of a part id, which selects a given B-Rep structure, and a selected entity id, which represents a piece of topology within the specified B-Rep structure.   Once a given topological entity is identified and retrieved from the B-Rep structure, one can navigate the connectivity information contained in the B-Rep to find other topological entities nearby.  By using a combination of the Clear() function and the AddToSelection(…) function associated with the selection manager, one can programmatically control which topological entities are currently selected. 

You can see this technique in action by perusing the Example in the previous posting.

An ANSYS Mechanical Scripting Example

Creating Unique Named Selections for Every Face in A Body.

I’ve put together a simple example of creating a unique named selection for each face associated with a selected set of bodies.  This example demonstrates the some of the techniques described in Part 3 of this series.  I recommend you read that post first to get a feel for the selection manager and ANSYS B-Rep data structures before returning to this post.

The unique name for each face is built up from a base name with the part ID and face ID appended to this name.  The function call used to create the named selection was taken directly from the ANSYS source code.  I used the technique described in Part 2 to find this function.  That is, I searched first for the string “Create Named Selection” and then for the string “ID_ComponentAdd” which led me to this XML code in dstoolbar.xml

   1:  <contextToolbar handle="NamedSelections">
   2:      <display id="2500">
   3:        <entry type="button" tipid="ID_ComponentAdd" icon="addselectiontolist">
   4:          <visibilityFilter name="idMeshing"/>
   5:          <visibilityFilter name="idMeshingAlbionSimAddin"/>
   6:          <visibilityFilter name="idAlbionSimAddin"/>
   7:          <actionCallback objectId="1" methodName="doAddComponent"/>
   8:          <visibilityCallback objectId="11" methodName="CanAddComponent"/>
   9:        </entry>

Which, as you can see on line 7, there is a method call named doAddComponent().  If we look into the call for this function, we see that it calls a function called addNamedSelection(bDialogFlag, specifiedName, criteria).  We can use this function to create a named selection with a programmatically generated name.

The flow of the example code is as follows:

  1. Check to see that at least one body is selected, if not, warn the user and quit.
  2. If one or more bodies are selected, build up a list of face ids and part ids associated with each body.
  3. Clear the selection manager so that nothing is selected.
  4. Iterate over the structure of part ids and face ids constructed in step 2 .
  5. For each face, create a unique name that is comprised of the part id and the face id concatenated to a base string.
  6. When finished, clear the selection manager so that nothing is selected.  This is sort of a clean up operation.  One enhancement might be to remember the bodies that were selected and clear out the faces from the selection manager and reselect these bodies so that the program is in the same state after this macro as it was when the macro was invoked.

Here is the code.  You can copy and paste this into a text file and save it somewhere on your hard drive.  The, select a body and choose Tools->Run Macro and navigate to the text file to see it work.

/****************************************************************************
 *                                                      Phoenix Analysis And
 *                                                       Design Technologies
 ****************************************************************************
 *  This script is designed to be run from the Tools->Run Macro menu from
 *  within ANSYS Mechanical.  The purpose of the script is to create a named
 *  selection for each individual face within the currently selected body.  
 *  This script makes use of the addNamedSelection() ANSYS Mechanical function 
 *  to actually create the component once the selection has been made.
 ***************************************************************************/

//  NOTE: The base name below is used as the first part of each named
//  selection name.  Change this to something more appropriate if you desire.
var basename = 'faces';

//  Here are a series of global variables that hold references to some
//  of the main objects in the mechanical application.

//  NOTE: The DS object is defined by the script engine ahead of time.  This
//  object references the main ANSYS Mechanical application

//  The script code object is a reference to the javascript engine embedded
//  inside the mechanical application.  This object has all of the functions
//  that are defined in the ANSYS mechanical javascript files
var SC = DS.Script;

//  The selection manager object is a reference to the ANSYS mechanical 
//  geometry selection tool.  This is a complex object that contains 
//  references to the currently selected set of geometric entities.  More
//  importantly, it contains back pointers to entity ids within the 
//  Mechanical BREP structure.
var SM = SC.sm;

//  NOTE:  The ANSYS Mechanical BREP structure is a classical topological
//  data structure that defines the connectivity of the geometry within
//  the currently opened model.  This structure consists of lists of
//  vertices, edges, faces, shells, bodies and parts.  More importantly,
//  the connectivity between these entities is represented in the data
//  structure as "use" lists.  For example, a vertex contains a list
//  of edges to which it is attached.  Likewise, an edge contains a list
//  of faces to which it is attached, as well as a pointer to the start
//  and end vertices of the edge.  With this information, one can "walk"
//  the topology of the part by iterating over the various use lists.
//  For this macro, we will grab all of the faces attached to a given
//  body using the BREP structure as a guide.

//  These type constants are encoded into the top two bits of a given
//  brep entity id.
var BODY_TOPO_TYPE = 3;
var FACE_TOPO_TYPE = 2;
var EDGE_TOPO_TYPE = 1;
var VERT_TOPO_TYPE = 0;

//  There is no concept of "main" in javascript, but it makes it easy to
//  find the entry point to the script.  There is a call to this function at
//  the bottom of this script
function main() {
    // First, see if we have a body selected.  If not, report to 
    // the user that a body must be selected first.
    if (SM.Parts > 0) {
        // Do the real work
        create_faces_ns_from_body();
    } else {
        SC.WBScript.Out('You must select a body before running the \
FacesFromBody.js macro.  Please select a body and re-run this \
macro.', true);
    }
}

// This function simply picks off the top two bits of a given topo id
// and returns the decimal value [0, 1, 2 or 3]  This corresponds to
// the topological types:
//      BODY_TOPO_TYPE = 3
//      FACE_TOPO_TYPE = 2
//      EDGE_TOPO_TYPE = 1
//      VERT_TOPO_TYPE = 0
function classify_entity_type(topo_id) {
    // The top two bits store the topological entity type within the topo
    // id value
    var type = topo_id >>> 30;
    return type;
}



// This function creates a named selection of all the faces associated with
// the currently selected body(s).  
// Pre Condition:   At least one body must be selected.
// Post Condition:  The selection manager will be cleared and a new named
//                  selection will exist for each face in the body.  The
//                  named selection will be of the form:
//                  basename_partid_faceid
function create_faces_ns_from_body() {
    // See structure definition below.
    var face_id_map = new Array(SM.Parts);

    // First we want to iterate over the selected parts and create
    // a list of all of the face ids for each selected part
    for (var i = 1; i <= SM.SelectedCount; i++) {
        var topo_id = SM.SelectedEntityTopoID(i);
        var part_id = SM.SelectedPartID(i);

        // Make sure this is a body.  If not, just continue around
        if (classify_entity_type(topo_id) != BODY_TOPO_TYPE) {
            continue;
        }
        var part = SM.PartMgr.PartById(part_id);
        var brep = part.BRep;
        // Pull the body object out of the BRep structure.  The call
        // for this is the Cell function with the topological id for
        // the body of interest
        body = brep.Cell(topo_id);
                
        // This array will be used to hold a list of face ids for a given
        // body.
        var face_ids = new Array();
        // These are the actual face objects associated with this body
        var faces = body.Faces;
        // Now store all of the face ids in our face id list
        for (var f = 1; f <= faces.Count; f++) {
            face_ids.push(faces(f).Id);
        }
        // Now that we have the face ids, put them into an array of structures
        // that each look like the following:
        // |-------|
        // |   0   |-> Part ID
        // |-------|
        // |-------|
        // |   1   |-> List of face ids for this part ID
        // |-------|
        face_id_map[i - 1] = new Array(2);
        face_id_map[i - 1][0] = part_id;
        face_id_map[i - 1][1] =
            face_ids.slice(0, face_ids.length); // Slice creates a copy of the array

    }
    // Now that we've built up our datastructure of face ids, we need to select them all
    SM.Clear();
    var name = null;
    for (var i = 0; i < face_id_map.length; i++) {
        var part_id = face_id_map[i][0];
        var face_ids = face_id_map[i][1];
        for (var j = 0; j < face_ids.length; j++) {
            SM.Clear();
            // Create a unique name based on the part id and face id
            name = basename + '_' + part_id.toString() + '_' + face_ids[j].toString();
            SM.AddToSelection(part_id, face_ids[j], false);
            // Create the component
            SC.addNamedSelection(false, name, SC.id_NS_UnknownMultiCriterion);           
        }
    }
    // Clear out the selection manager
    SM.Clear();
}

// Make it happen...
main();

ANSYS Mechanical Scripting: HOWTO Part 2

Where We Left Off Last Time…

In the first part of this how to guide, I left you with a general description of the architecture of ANSYS Mechanical product.  If you haven’t read that post, I recommend that you at least skim through it to familiarize yourself with the structure of ANSYS Mechanical.   One thing we established is that customizing ANSYS Mechanical is not for the faint of heart, but at the same time, it is also not an impossible task.  There are two primary difficulties associated with Mechanical customization.  They are:

  1. ANSYS Mechanical was never truly designed to be customized, at not least by end users.
  2. There is no documentation of the internal API, so we’re left with deciphering the inner workings ourselves.

Nonetheless, the Mechanical application is architected as a collection of functional pieces written in C++ that are tied together with Javascript via a COM interface.  What that means for us is that most of the functionality of ANSYS Mechanical is exposed to us through an interface that can be scripted.  So, even though Mechanical may not have been explicitly designed to be scripted nor are the internals documented, yet we find that it is possible to script, and since the ANSYS developers have followed a consistent pattern, its internal structure is comprehensible by an outside agent.

A Recipe for Discovering ANSYS Mechanical Functionality

As I hinted in the introduction to this post, the ANSYS developers have followed a sound software development practice in the construction of the ANSYS Mechanical product.  The practice is this:  Find a particular pattern that solves a problem you currently face in the most general case.  So, lets describe the problem the ANSYS developers faced and then we will describe what I think was their solution.   The problem can be summarized as this:  How do you create a GUI that is easily extended to later include functionality that you don’t even know about today?  That is the problem from 50,000 ft.  Here are some more details about the problem.

  • Extensible means that I (the ANSYS Developer) may need to add / subtract menu items, toolbar buttons, menu groups, button groups at will as I add in more functionality over time.  I don’t want to have to rewrite lots of low level code every time I need to add a button, for example.
  • I already know we (ANSYS) have customers around the world who speak different languages.  I don’t want to maintain 10 different versions of ANSYS Mechanical for every language that I need to support.  It would be nice if I could separate the language from the rest of the code.
  • If I’m going to be adding buttons, menu items, etc… at will, I need to be able to describe what happens when the user presses the button, selects the menu, etc…  I don’t want to paint myself into a corner, so this needs to be as general as possible.

So, that’s the problem and perhaps you’re already seeing the solution in your head.  Here is what the folks at ANSYS came up with, which is a pretty standard pattern for constructing a flexible GUI.

  1. Push the actual GUI construction to the very last minute.  That is, you don’t hard code a GUI, rather you hard code a GUI builder engine.  The GUI builder engine reads in a human readable text file at run time that describes the GUI layout and then it dynamically creates the GUI on the fly every time the program launches.   Now what you have is a static executable that can dynamically morph its appearance to anything you want.
  2. This is a collorary to 1.  Since we want to push GUI construction to the very last minute, we therefore have to have the means to describe at runtime what should happen when the user presses a button or selects a menu item in the most general terms.  Hence we need a dynamic scripting language that isn’t interpreted until run time.  Note that the only real purpose of this language is to orchestrate at a high level what should happen.  The low level details of say meshing for example can all be constructed in a compiled language and contained in libraries as long as we can call into those libraries using this scripting language.  So, because of 1 we end up with a scriptable interface to Mechanical.
  3. We want to isolate language differences, i.e. English, French, Germany, etc… to as small a subset of files as possible.  So, the best programming technique for this problem is the old “extra level of indirection.”  That is, instead of using a literal string like “File” in the code that describes the menu item “File”, we use an identifier ID_File that is just a numeric constant.  Then, we construct a table that maps a string value to a numeric identifier.  So, maybe ID_File has a numeric value of 438.  In our table there would be an entry of: 438, “File”.  Likewise, anywhere in the actual GUI that we need the string “File”, we substitute the constant ID_File in its place.  Now, lets say we want a French version of ANSYS Mechanical, what do we do?  We simply translate the string table entries into French and create a new string table file.  At run time, we read in the French string table and voilà, we have a French version of ANSYS Mechanical.

So, that is the approach the ANSYS developers took to solve the problem of creating an easily extensible GUI.  The benefit for us is that we get a scripting language to boot, and if we “run the engine backwards” we get an easy algorithm for figuring out what is going on under the hood in Mechanical.  So, here is the recipe for determining functionality within ANSYS Mechanical.

Main Ingredients

  1. Good Text Search Tool (I use Slickedit.  An added benefit is if it can search entire directories recursively)
  2. Patience

Procedure

  1. Locate the desired functionality in the GUI itself.
  2. Search the file dsstringtable.xml for the menu string, or toolbar button description string
  3. Note the corresponding ID_*** string ID within the string table file.
  4. Search for the ID_*** in dscontextmenu.xml or similar GUI structure file.
  5. When you find the ID_***, note the do*** command listed in the section associated with that GUI element’s action callback.
  6. Recursively search the aisol, or DesignSpace directory for the corresponding do*** command.

Could it possibly be simpler?

Example Use of the Above Recipe

I want to add a command object to a particular body in the tree, and eventually I want to do it programmatically.

Step 1: Locate the Functionality in the GUI Itself

Below is a screen shot of the context menu hierarchy one navigates through to insert a command object on a body in the tree.  So, we see that the string we are looking for is “Commands”

image

Step 2: Search dsstringtable.xml for the Menu String: “Commands”

Note that the string tables are located underneath a language directory structure.  If you haven’t currently guessed based on the content of this blog post, we’re using the string table located under Language\en-us\xml.  When we perform the search we find the following:

image

Step 3: Note the Corresponding ID_*** in the String Table

You can see from the image above that the ID associated with the string Commands is ID_CommandEditorDefaultName.  Sounds plausible.

Step 4: Search for ID_CommandEditorDefaultName in dscontextmenu.xml

Note, that the structure of the ANSYS Mechanical GUI is described in a series of xml files as mentioned in the sections above.  The file dscontextmenu.xml is one such file that shockingly describes the layout of the various context menus within the program.  So, lets search for the string ID_CommandEditorDefaultName in that file and see what we find.  Here are the results of that search:

image

Note the structure of this file.  There are various entries that represent individual menu items within a given context menu.  For each entry there are a series of action callbacks that have an associated methodName attribute.  By deduction we can assume that when the user clicks on this menu item, the listed javascript functions are called to perform the menu item’s action.  You’ll note too that there is a visibilityCallback as well.  Those wacky ANSYS developers thought of everything…

Now, one of the action callbacks is called “doPrototypeInsertCommandEditor”.  If that doesn’t sound like a function that will insert a command editor on a prototype, I don’t know what does.  Let’s see if we can find it.

Step 5: Search for doPrototypeInsertCommandEditor in the aisol or Designspace Directories

Now we search inside the javascript (.js) files for a function called doPrototypeInsertCommandEditor() and see what we can find.  Inside the file DSMenuScript.js we find this function definition:

image

Well, what do you know?  There is the 21 line function that gets called whenever you insert a command object onto a body in tree.   In the intervening time between this post and the next, spend some time staring at this code.  There are some fundamental aspects of the architecture of ANSYS Mechanical that appear even in this short function.  Now, I know your probably thinking to yourself, this is all fine and dandy, but what if I don’t just want to insert a command object.  What if I want the command object to have some APDL in it?  Furthermore, what do all those other variables, etc… really mean.  Lastly, how do I even call this blasted function?  Stay tuned… there is much, much more to come.

ANSYS Mechanical Scripting: HOWTO Part 1

Background

Over the past 10+ years, ANSYS has steadily migrated more and more functionality into the DesignSpace, wait Workbench… no I mean Workbench Simulation or rather now ANSYS Mechanical product line.  If you’ve used ANSYS in the last 5 years you’ve probably come to know this interface as Workbench, and you’ve probably delineated it in your mind as standing in opposition to the old interface, “Classic” or Mechanical APDL.  Perhaps these qualities of each interface come to mind for you as they do for me

Mechanical APDL is:

  1. Old looking and visually clunky.
  2. Quirky, but powerful.
  3. The only thing “real” analyst use.
  4. Scriptable, Scriptable and More Scriptable. 

ANSYS Mechanical is:

  1. Cool, slick looking user interface
  2. Super easy to use when your analysis is supported out-of-the-box.
  3. Great meshing and geometry support.
  4. Not scriptable.

Fortunately for us the margin in capabilities between Mechanical APDL and ANSYS Mechanical is narrowing at each release and so more and more “real” analyst are turning to ANSYS Mechanical for their simulation needs.   However, one glaring difference between the two products that has remained virtually unchanged since the beginning is in the area of scripting.   Mechanical APDL is scriptable at its core.  In fact, that is the only way I personally interact with the program, which I’m sure is true for many of you as well.  ANSYS Mechanical on the other hand appears as though it is just an impenetrable black box.  It does what it does and that is that.  Fortunately, that is only partially true.  The reality is there is a whole underworld that exists in ANSYS Mechanical that is scriptable, it is just very cleverly and discretely hidden.  (i.e. It isn’t documented)  ANSYS, Inc. plans to keep it that way the best I can tell.  However, I’ve spent the last few months of my life painstakingly pulling back the layers and peeking inside.  So, I hope to show you some interesting tidbits I’ve found along the way.

Anatomy of DesignSpace

Why am I going back to calling ANSYS Mechanical DesignSpace?  Well, as you might imagine, developers aren’t quite as flexible as the marketing types, so a lot of the code that is ANSYS Mechanical still carries the moniker DesignSpace, or DS.  You can imagine the mutiny that would occur, if every time a suit got the bright idea to change a product name, all the developers were required to rewrite their code to reflect the new name.  It ain’t happening.

So, what does DS look like on the inside?  Well, the best I can tell it is architected this way:

ApplicationArchitecture

The outer shell (aka the impenetrable black box) is the GUI.  Obviously, this is how the user typically interacts with the program.  However, the GUI itself is not a statically compiled piece of executable code.  It is more like an interpreter that builds itself every time you launch ANSYS Mechanical. (You may be asking yourself, “Self, is this why it takes so bloody long to load…”  Self replies “Perhaps…  BTW, I need more coffee.”)  So, the next logical question is “How does it know what to build?”  I’m glad you asked, because this turns out to be one of the keys to unlocking the black box.  The structure of the ANSYS Mechanical GUI is described in a handful of XML files.  This is important for a couple of reasons:

  1. XML is just text, so you don’t have to put on your binary glasses to read it.
  2. XML has structure. Brains really like structure.

It may seem like I’m getting nowhere fast with this, but hang with me.  Where we are at now is that the structure of the GUI is described by an XML file.  So what that means is that the text associated with all of the menus and buttons in the GUI; whether a GUI entity is a toolbar, or a toolbar button or a menu item; all of that is described textually inside these XML files.  However, a GUI that doesn’t interact with the user is just a G.  How does the UI come into play? 

I’m glad you asked that as well.  The interactive part of the ANSYS Mechanical UI is handled primarily through the use of javascript.  Javascript is an interpreted scripting language that is most popular in the context of web development.  It is text based, i.e. it isn’t compiled.  The way it works is that there is a javascript interpeter built into the ANSYS Mechanical executable.  Almost all of the core functionality of ANSYS Mechanical is implemented as a set of C++ libraries.  These are just great big pieces of code like meshing, or virtual topology, or graphics or geometry selection that expose a set of routines that hook into the javascript interpreter.  All of these functional pieces are then glued together with javascript.  Thus, it turns out that the UI part of the ANSYS Mechanical GUI is just a whole bunch of javascript code juggling all of these big pieces of functionality that are implemented in C++.  It is not completely unlike Mechanical APDL, where a function like ASBW hides a complex surface-to-surface intersection routine followed by a BREP patchup; however, we simply know it as “Divide these areas by the working plane.”   In the same way there are various javascript calls that a script can issue that will insert a new boundary condition, for example. Or, perhaps, change a mesh control and re mesh.

Finally, I’ll leave you with one last piece of the puzzle for this post, then we’ll pick it up with examples in the next post.   You remember that I said that the GUI is cooked up each time using an XML file as the recipe.  Inside that XML file are definitions for javascript functions that are called whenever a user clicks on a button, or selects a menu item.  This provides the link between the static GUI structure and the dynamic user interaction.  This is how we turn the key to unlock the black box.  By searching in the xml file to find a GUI object that does what we want, we can then determine the javascript function that is called when the user interacts with that GUI object.  By finding the javascript function name, we can find the function implementation inside the ANSYS installation directory.  By finding the implementation, we can study it and figure out what it does. 

Stay tuned for the next installment to see this in action.