How to Use the ModelPart

Revision as of 14:30, 26 November 2007

The model part is the place were all the informations about one part of the model are stored.

It is of course possible to access directly to the "arrays" of Nodes, Elements, Conditions and so on. This can be done by typing

  ModelPart::NodesContainerType& rNodes = model_part.Nodes(); //gives access to the list of nodes
  ModelPart::ElementsContainerType& rElements = model_part.Elements(); //provides access to all of the Elements
  ModelPart::ConditionsContainerType& rConditions = model_part.Conditions(); //provides access to all Conditions

All of these containers are designed to allow "stl style" operations. The implemenation is however non standard and an effective use requires some understanding of their design concepts.

Design Concepts

The basic problem is that in the finite Element context a non-contiguous numeration of the Nodes (or Elements) is often encountered, and has to be dealt with. In the Kratos, the NodesContainerType implements an ordered vector in which the entries are stored following their "Id".

To make an example, let's consider a model part containing 4 nodes with the following IDs and coordinates

It is immediately evident that some IDs are missing from the list, making unconvenient for example to store the node "10" at the position "10" of an array. In the Kratos the choice is made to store exclusively the Nodes that effectively exist, which in practice implies that the four nodes will be stored at the positions

 0 1 2 3

of the "NodesContainer".

As a design concept the data structure used for storing the objects is somehow hidden to the final user and is subjected to possible future changes (this is the reason to introduce the typedef "ModelPart::NodesContainerType" ). Some operations will be however guaranteed to be efficient while others will be kept "as efficient as possible" without being optimal.

Iterating over all the objects in the container is designed to be the most efficient operation. As an example the loop

 for( 
     ModelPart::NodesContainerType::iterator inode = model_part.NodesBegin();
     inode != model_part.NodesEnd();
     inode++)	
 {
     KRATOS_WATCH(inode->Id()); //print the value of the node Id
 }

print the Id of all the nodes ... in the most efficient way.

Operators [] , () and find()

At this point a question naturally rises: Can i access a node by Id? The answer is YES and in a resonably efficient way. The operator [] is overloaded so to give access by Id, this implies that the code

  unsigned int X = model_part.Nodes()[10].X(); //get the id of Node 1
  KRATOS_WATCH(X); // prints the  coordinate X of the node with Id 10

provides the output "12.33" even if the node is stored (following the example described above) in the 4th position of the list This notation is consistent with the definition of the "set" in the standard library, which allows to apply directly all of the stl algorithms.

We should however be aware that if the node Id does not exist in the list (Implying that the node is not a part of the ModelPart), any attempt to accessing to the node will lead to a runtime error.

The operator [] returns a reference to the node with the corresponding Id. The operator () has an identical behaviour but returns a Pointer to the same node. As an example

  Node<3>& rnode = model_part.Nodes()[10]; //this provides a reference to the node with Id 10
  Node<3>::Pointer = model_part.Nodes()(10); //this one provides a Pointer to node 10

The possibility exists to look for a node by index without reasking a runtime error for non-existing nodes. This is achieved by the operator "find"

  ModelPart::NodesContainerType::iterator inode = model_part.Nodes().find(3) 

which returns an iterator pointing to the node we were looking for, or, if the node was not found, pointing to the "model_part.Nodes().end()"

Elements and Conditions

The example shown above concentrates on Nodes. The interface for "Elements" and "Conditions" is very similar. Iteration by index is identical (with the same problems). Iteration by iterators is

 for( 
     ModelPart::ElementsContainerType::iterator iel = model_part.ElementsBegin();
     iel != model_part.ElementEnd();
     iel++
     )	

for elements, and

 for( 
     ModelPart::ConditionsContainerType::iterator icond = model_part.ConditionsBegin();
     icond != model_part.ConditionsEnd();
     icond
     )	

for conditions

Alternatives for Iteration

As a matter of fact the NodesContainerType is currently implemented as a Vector (guaranteed to be ordered following the Ids). In some specific cases it may be convenient to access directly to the node by knowing its position in the underlying container rather than by Id.

For example the Id of all of the nodes can be plotted using this alternative approach

  for(unsigned int i = 0; i!=model_art.Nodes().size(); i++)
      KRATOS_WATCH((model_part.NodesBegin()+i)->Id());

in which an integer "i" is used as counter and the access is obtained by random access iterator.

Apart for this the possibility exist to iterate directly on the underlying container (accessible by "model_part.Nodes()).GetContainer()"). This choice should not be performed as it depends on the actual implementation of the container and may lead to errors if the behaviour of the NodesContainer is changed in the future.

Filling a new model_part

An operation that is often needed is to generate a new model part given another one, for example to represent a different physical behaviour. As an example one may want to take the mesh used to represent a structural problem and to generate a set of new elements with the same connectivity but representing a thermal problem.

In the practice it is common to have two (or more) model parts sharing the same nodes (and consequently the same database) but having different elements to discretize the continua.

Let's start considering that a new_model_part is created and that we want to use the same nodes as for the "existing_model_part". The "copy" of the array of nodes can be done by simply typing

  new_model_part.Nodes() = existing_model_part.Nodes();

this ensures that the new model part is filled with pointers to the nodes that were before in the existing_model_part. Note that no copy of the nodes is made and the new model part still uses the same nodes (now shared between the two models). The command simply implies filling an array of pointers.

In other cases it may be convenient to fill the array of nodes of the "new_model_part" only with a part of the nodes of the existing one. The following code fills the new model with all of the nodes having Id bigger than 100

  for(NodesContainerType::iterator it = existing_model_part.NodesBegin(); it!=existing_model_part.NodesEnd(); it++)
  {
     if(it->Id() > 100)
        new_model_part.push_back( *(it.base() ) ); //take care!! this is how to get a pointer having the iterator!
  }

We want now to fill the new model part with elements of type "MyNewElementType" by reusing the connectivity of the existing model part.

  for(ElementsContainerType::iterator it = existing_model_part.ElementsBegin(); it!=existing_model_part.ElementsEnd(); it++)
  {
      Geometry< Node<3> >::Pointer pgeom = it->pGetGeometry(); //get the connectivity of the origin element
      unsigned int new_id = it->Id(); //we want to create a new element with the same Id as the old one
      PropertiesType::Pointer pprop = it->GetProperties(); //get the properties (of the old model part)
      //create the new element by using the old nodes the same Id as the origin element and the same properties (owned by the old model part)
      new_model_part.Elements().push_back( Element::Pointer( new MyNewElementType(new_id,pgeom,pprop) ) );
  }

as an alternative we can fill all of the new model part by using a given type of elements we know by name. Say for example "BaseTemperatureElement". The idea in this case is to access to a "model element" which the kratos has stored, and clone it at will.

  //get the "model element" named "BaseTemperatureElement"
  Element const& r_reference_element = KratosComponents<Element>::Get("BaseTemperatureElement"); 
  //loop to make the copy
  for(ElementsContainerType::iterator it = existing_model_part.ElementsBegin(); it!=existing_model_part.ElementsEnd(); it++)
  {
      Geometry< Node<3> >::Pointer pgeom = it->pGetGeometry(); //get the connectivity of the origin element
      unsigned int new_id = it->Id(); //we want to create a new element with the same Id as the old one
      PropertiesType::Pointer pprop = it->GetProperties(); //get the properties (of the old model part)
      //create a copy using the "reference element" as a model
      Element::Pointer p_element = rReferenceElement.Create(new_id, pgeom ,properties);
      //fill the new model part 
      new_model_part.Elements().push_back(p_element);
  }

Common Programming Errors

All of the implementation used in the Kratos Container relies rather heavily on the concepts of the C++ standard library. One of the important points to be understood is the difference between the actual "size" of a container and its "capacity" which is the memory actually allocated for it. We do not want to discuss this aspect here in any detail, nevertheless the user should be aware that once one of the model_part's containers is allocated more memory is reserved than the one that is actually used and that this memory is used to create a "buffer" with the aim of speeding up operations as the push_back (adding an element to the end). The idea is basically that one will be able to push_back efficiently to a vector as long as its capacity is not exceeded, while at the moment at which this happens the vector is reallocated somewhere else in the memory and the capacity is increased (typically doubled at every reallocation). This reallocation may lead to unexpected problems which we will try to explain with a simple example:

let's assume we have an model_part with 1000 elements and we want to add them one by one to another model part iterating "by integer". Let's consider the code

  ElementsContainerType::iterator destination_begin = destination_model_part.ElementsBegin();
  ElementsContainerType::iterator origin_begin = origin_model_part.ElementsBegin();
  for(unsigned int i = 0; i!=1000; i++)
  {
     //save the elements of the origin model part as members (as well) of the new model part
     destination_model_part.Elements().push_back( *( (origin_begin + i).base() );
     //
     //change the Id (note that this will change both the origin and destination elements as we are using pointers
     (destination_begin + i)->Id() = i;  //POTENTIAL SEGMENTATION FAULT!!!!!!
  }

this example will most likely lead to run time errors because of the line

     (destination_begin + i)->Id() = i;  

the problem is here that the container of the elements of the destination model part is most likely reallocated sometimes while pushing back the 1000 pointers. When this happens the reference to the beginning of the array will be invalidated as the array will be copied somewhere else in the memory, and the former position will be freed. An attempt to access to "(destination_begin + i)" will thus "break" the memory and will be invalid. The cure for this is very simple as

     //(destination_begin + i)->Id() = i;  //WRONG!!

should be replaced to

     //(destination_model_part.ElementsBegin() + i)->Id() = i; //OK!! 

the second one is safe at the "ElementsBegin" is automaticall updated when the vector is reallocated