Personal tools
You are here: Home / cm / Wiki / Write out material points positions
Navigation
Log in


Forgot your password?
 

Write out material points positions

This is a tip in case you would like to know the coordinates of material
points in your mesh.

My motivation was that I wanted to track the position in space of many material points in my mesh as it deformed through time. The solution for my problem was to use the command 'fem evaluate field'. (Other options may be to use commands for listing stresses or using grid points which I haven't tried out but I assume the coordinates may be printed out for these commands too.)

The fem evaluate field command listed out xi-positions and rectangular Cartesian coordinates and a field value. (I used RC-coordinates. It may print out other types of coordinates if other coordinate system is used.) The limitation of this command is that there will be some duplication of points (on element interfaces) and it will not produce a regular grid in RC coordinates.

Below follows part of my script that first prints out the field-file, then a subroutine extracts the coordinates from this file and writes them to a new file (together with a parameter that says if the point is on the element boundary):

$xi_res=0.2;  # chosen resolution

$xi1_res=$xi_res;

$xi2_res=$xi_res;

$xi3_res=$xi_res;

$FILE = "pointFiles/tempList";

fem evaluate field;$FILE xi $xi1_res,$xi2_res,$xi3_res reduced dependent;

$return = MyFunc::writePts("$FILE.opfiel",$newTime);





# Subroutine in MyFunc.pm:

sub writePts

{

  my ($inFile, $time) = @_;



  my $dir = "pointFiles/";

  my $outFile = sprintf("./$dir/points_time_%04d.txt", $time);



  CORE::open(INFID, $inFile) || die "Error: Couldn't open $inFile\n";

  CORE::open(OUTFID, "> $outFile") || die "Error: Couldn't open $outFile\n";

  while (my $line = <INFID>)

  {

    if($line =~ m/ Xi:  (.+) (.+) (.+)  Z:\s+(.+)\s+(.+)\s+(.+)\s+(.+)/)

    {

      my $bPrintPt = 1;

      my $xi1 = $1;

      my $xi2 = $2;

      my $xi3 = $3;

      my $x = $4;

      my $y = $5;

      my $z = $6;

      my $bSurf = 0;



      #Check if point may be a duplicate point, i.e. is a surface point

      if(($xi1==0) || ($xi2==0) || ($xi3==0) || ($xi1==1) || ($xi2==1) || ($xi3==1))

      {

        $bSurf = 1;

      }

      if($bPrintPt)

      {

        print OUTFID sprintf("%d  %9.4f  %9.4f  %9.4f\n", $bSurf, $x, $y, $z);

      }

    }

  }

  CORE::close(INFID);

  CORE::close(OUTFID);

  return(1);

}

I also wrote a script to remove duplicated points from the coordinate files. The removal process is very slow if there are many points in each element and could very likely be optimized. As I needed this done for 100 files, I generated two routines: the first removed the duplicated points in file 1 and kept track of the line numbers with non-duplicated points. In file 2 to 100 it only read in those lines, which sped up the process dramatically for those files. See scripts below:

$inFile = "pointFiles/points_time_0001.txt";

$return = writePts1($inFile, 1);





for($time=2;$time<=100;$time++)

{

  $inFile = sprintf("pointFiles/points_time_%04d.txt", $time);

  print $inFile, "\n";

  $return = writePts2($inFile, $time);

}







sub writePts1

{

  my ($inFile, $time) = @_;



  my $dir = "pointFiles/";

  my $outFile = sprintf("./$dir/materialPoints_time_%04d.txt", $time);



  my $maxDupPts = 100000; # An initial guess to allocate memory

  my $counter=0;

  my %dupPts_x; #array of points that may be duplicated, i.e. xi value of 0 or 1

  my %dupPts_y; #array of points that may be duplicated, i.e. xi value of 0 or 1

  my %dupPts_z; #array of points that may be duplicated, i.e. xi value of 0 or 1

  $#dupPts_x = $maxDupPts;

  $#dupPts_y = $maxDupPts;

  $#dupPts_z = $maxDupPts;

  my $dupCounter = 0;

  my $lineNumber = 0;

  my $lineFile = "./$dir/writtenLines.txt";

  CORE::open(NUMFID, "> $lineFile") || die "Error: Couldn't open $lineFile\n";

  CORE::open(INFID, $inFile) || die "Error: Couldn't open $inFile\n";

  CORE::open(OUTFID, "> $outFile") || die "Error: Couldn't open $outFile\n";

  while (my $line = <INFID>)

  {

    $lineNumber++;

    if($line =~ m/([0,1]) \s+ (.+) \s+ (.+) \s+ (.+)/)

    {

      my $bPrintPt = 1;

      my $bSurf = $1;

      my $x = $2;my $outFile = sprintf("./$dir/junk_time_%04d.txt", $time);

      my $y = $3;

      my $z = $4;

      #print sprintf("b %d  x%9.4f  y%9.4f  z%9.4f\n", $bSurf, $x, $y, $z);

      #Check if point is a duplicate

      if($bSurf)

      {

        for(my $n=0;$n<$dupCounter;$n++)

        {

          if(($x==$dupPts_x[$n]) && ($y==$dupPts_y[$n]) && ($z==$dupPts_z[$n])) #Point has been printed

          {

            $bPrintPt=0;

            last;

          }

        }

        if($bPrintPt) #Point has not been printed: add to list

        {

          $dupPts_x[$dupCounter] = $x;

          $dupPts_y[$dupCounter] = $y;

          $dupPts_z[$dupCounter] = $z;

          $dupCounter++;

          print $dupCounter, "\n";

        }

      }

      if($bPrintPt)

      {

        print OUTFID sprintf("%9.4f  %9.4f  %9.4f\n", $x, $y, $z);

        #print sprintf("%14.9f  %14.9f  %14.9f\n", $4, $5, $6);

        $counter++;

        print NUMFID "$lineNumber\n";

      }

    }

  }

  if($dupCounter > $maxDupPts)

  {

    print "\nWARNING: maxDupPts variable should be increased!\n";

  }

  CORE::close(INFID);

  CORE::close(OUTFID);

  CORE::close(NUMFID);

  return(1);

}





sub writePts2

{

  my ($inFile, $time) = @_;



  my $dir = "pointFiles/";

  my $outFile = sprintf("./$dir/materialPoints_time_%04d.txt", $time);



  my $lineNumber = 0;

  my $lineWithNextPt = 0;

  my $lineFile = "./$dir/writtenLines.txt";

  my $bMorePts = 1;

  CORE::open(NUMFID, "$lineFile") || die "Error: Couldn't open $lineFile\n";

  CORE::open(INFID, $inFile) || die "Error: Couldn't open $inFile\n";

  CORE::open(OUTFID, "> $outFile") || die "Error: Couldn't open $outFile\n";

  unless($lineWithNextPt = <NUMFID>)

  {

    die("Error: Did not find file or line number in file $lineFile!\n");

  }

  while((my $line = <INFID>) && ($bMorePts))

  {

    $lineNumber++;

    #print "Next line to read: ", int($lineWithNextPt), "\n";

    #print "$lineNumber\n";

    if($lineNumber == int($lineWithNextPt))

    {

      if($line =~ m/([0,1]) \s+ (.+) \s+ (.+) \s+ (.+)/)

      {

        my $bSurf = $1;

        my $x = $2;

        my $y = $3;

        my $z = $4;

        print OUTFID sprintf("%9.4f  %9.4f  %9.4f\n", $x, $y, $z);

        #print sprintf("%14.9f  %14.9f  %14.9f\n", $4, $5, $6);

      }

      unless($lineWithNextPt = <NUMFID>)

      {

        $bMorePts = 0;

        #$lineWithNextPt = int($lineWithNextPt);

        #print "Next line to read: ", int($lineWithNextPt), "\n";

      }

    }

  }

  CORE::close(INFID);

  CORE::close(OUTFID);

  CORE::close(NUMFID);

  return(1);

}