#include "roccom.h"
#include <cstdio>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <cstring>
#include <string>
#include <cstdlib>
#include <cmath>
#include <cassert>
#include <sstream>
#include "roccom_assertion.h"
#include "../Rocsurf/test/IM_Reader.h"

Include dependency graph for rotate.C:

Classes
struct	Control_parameter

Functions
	COM_EXTERN_MODULE (Rocblas)

	COM_EXTERN_MODULE (Rocmap)

	COM_EXTERN_MODULE (Rocprop)

	COM_EXTERN_MODULE (Rocsurf)

	COM_EXTERN_MODULE (Rocout)

void	load_modules ()

void	print_usage (int argc, char *argv[])

void	read_control_file (const char *fname, Control_parameter &cp)

void	init_parameters (const Control_parameter &cntr_param)

std::string	read_in_mesh (const char *fname)

void	init_constraints_acmfrac (const string &wname)

void	init_constraints_acmflu (const string &wname)

void	init_constraints_starslice (const string &wname)

void	init_attributes (const string &wname, const Control_parameter &cntr_param)

double	square (double x)

void	init_speed_arrays (const double coors, double velocities, int n)

void	init_speed (const string &wname)

void	output_solution (const string &wname, const char *timelevel)

int	main (int argc, char *argv[])

Variables
static int	rank = 0

Function Documentation

COM_EXTERN_MODULE ( Rocblas )

COM_EXTERN_MODULE ( Rocmap )

COM_EXTERN_MODULE ( Rocprop )

COM_EXTERN_MODULE ( Rocsurf )

COM_EXTERN_MODULE ( Rocout )

void init_attributes	(	const string &	wname,
		const Control_parameter &	cntr_param
	)

Definition at line 332 of file rotate.C.

References COM_CHAR, COM_DOUBLE, COM_INT, COM_new_attribute(), COM_resize_array(), COM_set_size(), and COM_window_init_done().

                                                            {
   COM_new_attribute((wname+".flag").c_str(), 'p', COM_DOUBLE, 1, "");
   COM_set_size( (wname+".flag").c_str(), 0, 1);
 
   COM_new_attribute((wname+".disps_total").c_str(), 'n', COM_DOUBLE, 3, "");
   COM_new_attribute((wname+".cnstr_types").c_str(), 'n', COM_INT, 1, "");
   
   COM_new_attribute((wname+".spds").c_str(), 'n', COM_DOUBLE, 3, "m/s");
   COM_new_attribute((wname+".disps").c_str(), 'n', COM_DOUBLE, 3, "m");
 
   COM_new_attribute((wname+".disps_novis").c_str(), 'n', COM_DOUBLE, 3, "");
   COM_new_attribute((wname+".facenormals").c_str(), 'e', COM_DOUBLE, 3, "");
   COM_new_attribute((wname+".faceheights").c_str(), 'e', COM_DOUBLE, 1, "");
 
   // Attribute for storing the number of eigenvalues for each node.
   COM_new_attribute((wname+".lambdas").c_str(), 'n', COM_DOUBLE, 3, "");
   COM_new_attribute((wname+".eigvecs").c_str(), 'n', COM_DOUBLE, 9, "");
   COM_new_attribute((wname+".tangranks").c_str(), 'n', COM_CHAR, 1, "");
   COM_new_attribute((wname+".scales").c_str(), 'n', COM_DOUBLE, 1, "");
 
   COM_resize_array( (wname+".atts").c_str());
   COM_window_init_done( wname.c_str());
 }

Here is the call graph for this function:

void init_constraints_acmflu ( const string & wname )

Definition at line 260 of file rotate.C.

References COM_call_function(), COM_get_array(), COM_get_attribute_handle(), COM_get_function_handle(), COM_get_size(), j, and x.

                                                    {
   static int BLAS_copy_scalar=0, cnstr_handle, MAP_maxabs;
   
   if ( BLAS_copy_scalar==0) {
     COM_get_function_handle( "BLAS.copy_scalar");
     cnstr_handle = COM_get_attribute_handle( (wname+".cnstr_types").c_str());
     MAP_maxabs = COM_get_function_handle( "MAP.reduce_maxabs_on_shared_nodes");
   }
 
   int zero = 0;
   COM_call_function( BLAS_copy_scalar, &zero, &cnstr_handle);
 
   // Set for pane 3 separately
   int *cnstr_types;
   COM_get_array( (wname+".cnstr_types").c_str(), 3, (void**)&cnstr_types);
   
   MAP::Vector_3<double> *coors;
   COM_get_array( (wname+".nc").c_str(), 3, (void**)&coors);
 
   int nnodes;
   COM_get_size( (wname+".nc").c_str(), 3, &nnodes);
 
   // Loop through all the points
   for ( int j=0; j<nnodes; ++j) {
     if ( coors[j].x() <= 1.e-10 || coors[j].x() > 0.047624)
       cnstr_types[j] = 2;       // Fix the nodes at two ends.
     else 
       cnstr_types[j] = 'x'; // Allow the nodes to move only along x dir.
   }
 
   COM_call_function( MAP_maxabs, &cnstr_handle);
 
   int PROP_set_cnstr = COM_get_function_handle( "PROP.set_constraints");
   COM_call_function( PROP_set_cnstr, &cnstr_handle);
 }

Here is the call graph for this function:

void init_constraints_acmfrac ( const string & wname )

Definition at line 226 of file rotate.C.

References COM_call_function(), COM_get_array(), COM_get_attribute_handle(), COM_get_function_handle(), COM_get_size(), and j.

                                                     {
   static int BLAS_copy_scalar=0, cnstr_handle, MAP_maxabs;
   
   if ( BLAS_copy_scalar==0) {
     COM_get_function_handle( "BLAS.copy_scalar");
     cnstr_handle = COM_get_attribute_handle( (wname+".cnstr_types").c_str());
     MAP_maxabs = COM_get_function_handle( "MAP.reduce_maxabs_on_shared_nodes");
   }
 
   int zero = 0;
   COM_call_function( BLAS_copy_scalar, &zero, &cnstr_handle);
 
   // Set for pane 2 separately
   int *cnstr_types;
   COM_get_array( (wname+".cnstr_types").c_str(), 2, (void**)&cnstr_types);
   
   MAP::Vector_3<double> *coors;
   COM_get_array( (wname+".nc").c_str(), 2, (void**)&coors);
 
   int nnodes;
   COM_get_size( (wname+".nc").c_str(), 2, &nnodes);
 
   // Loop through all the points
   for ( int j=0; j<nnodes; ++j) {
     cnstr_types[j] = 'x'; // Allow the nodes to move only along x dir.
   }
 
   COM_call_function( MAP_maxabs, &cnstr_handle);
 
   int PROP_set_cnstr = COM_get_function_handle( "PROP.set_constraints");
   COM_call_function( PROP_set_cnstr, &cnstr_handle);
 }

Here is the call graph for this function:

void init_constraints_starslice ( const string & wname )

Definition at line 297 of file rotate.C.

References COM_call_function(), COM_get_array(), COM_get_attribute_handle(), COM_get_function_handle(), COM_get_size(), j, x, and y.

                                                       {
   static int BLAS_copy_scalar=0, cnstr_handle, MAP_maxabs;
   
   if ( BLAS_copy_scalar==0) {
     COM_get_function_handle( "BLAS.copy_scalar");
     cnstr_handle = COM_get_attribute_handle( (wname+".cnstr_types").c_str());
     MAP_maxabs = COM_get_function_handle( "MAP.reduce_maxabs_on_shared_nodes");
   }
   int zero = 0;
   COM_call_function( BLAS_copy_scalar, &zero, &cnstr_handle);
 
   // Set for pane 2 separately
   int *cnstr_types;
   COM_get_array( (wname+".cnstr_types").c_str(), 2, (void**)&cnstr_types);
   
   MAP::Vector_3<double> *coors;
   COM_get_array( (wname+".nc").c_str(), 2, (void**)&coors);
 
   int nnodes;
   COM_get_size( (wname+".nc").c_str(), 2, &nnodes);
 
   // Loop through all the points
   for ( int j=0; j<nnodes; ++j) {
     if ( coors[j].y() > -20.054 || coors[j].x() < -20.5538)
       cnstr_types[j] = -'y';  // Move within xz plane.
     else
       cnstr_types[j] = 2;     // Fix other points
   }
 
   COM_call_function( MAP_maxabs, &cnstr_handle);
 
   int PROP_set_cnstr = COM_get_function_handle( "PROP.set_constraints");
   COM_call_function( PROP_set_cnstr, &cnstr_handle);
 }

Here is the call graph for this function:

void init_parameters ( const Control_parameter & cntr_param )

Definition at line 155 of file rotate.C.

References COM_call_function(), COM_get_function_handle(), Control_parameter::courant, Control_parameter::eigthres, Control_parameter::fangle, Control_parameter::method, Control_parameter::normaldif, rank, Control_parameter::rediter, Control_parameter::smoother, Control_parameter::verbose, and Control_parameter::wavefrontal.

                                                            {
   int PROP_set_option = COM_get_function_handle( "PROP.set_option");
 
   if ( !cntr_param.method.empty()) {
     COM_call_function( PROP_set_option, "method", cntr_param.method.c_str());
 
     if ( rank==0) std::cout << "Set propagation method to " << cntr_param.method << std::endl;
   }
 
   if ( !cntr_param.wavefrontal.empty()) {
     COM_call_function( PROP_set_option, "wavefrontal", cntr_param.wavefrontal.c_str());
     if ( rank==0) std::cout << "Set wavefrontal to " << cntr_param.wavefrontal << std::endl;
   }
 
   if ( !cntr_param.normaldif.empty()) {
     COM_call_function( PROP_set_option, "normaldif", cntr_param.normaldif.c_str());
     if ( rank==0) std::cout << "Set normaldif to " << cntr_param.normaldif << std::endl;
   }
 
   if ( !cntr_param.eigthres.empty()) {
     COM_call_function( PROP_set_option, "eigthres", cntr_param.eigthres.c_str());
     if ( rank==0) std::cout << "Set eigthres to " << cntr_param.eigthres << std::endl;
   }
 
   if ( !cntr_param.courant.empty()) {
     COM_call_function( PROP_set_option, "courant", cntr_param.courant.c_str());
     if ( rank==0) std::cout << "Set courant constant to " << cntr_param.courant << std::endl;
   }
 
   if ( !cntr_param.fangle.empty()) {
     COM_call_function( PROP_set_option, "fangle", cntr_param.fangle.c_str());
     if ( rank==0) std::cout << "Set feature angle to " << cntr_param.fangle << std::endl;
   }
 
   if ( !cntr_param.smoother.empty()) {
     COM_call_function( PROP_set_option, "smoother", cntr_param.smoother.c_str());
     if ( rank==0) std::cout << "Set smoother to " << cntr_param.smoother << std::endl;
   }
   
   if ( !cntr_param.rediter.empty()) {
     COM_call_function( PROP_set_option, "rediter", cntr_param.rediter.c_str());
     if ( rank==0) std::cout << "Set rediter to " << cntr_param.rediter << std::endl;
   }
 
   if ( !cntr_param.verbose.empty()) {
     COM_call_function( PROP_set_option, "verbose", cntr_param.verbose.c_str());
 
     if ( rank==0) std::cout << "Set verbose level to " << cntr_param.verbose << std::endl;
   }
 
 }

Here is the call graph for this function:

void init_speed ( const string & wname )

Definition at line 377 of file rotate.C.

References COM_free_buffer(), COM_get_array(), COM_get_panes(), COM_get_size(), i, init_speed_arrays(), and n.

Referenced by main().

                                       {
   // Set speed for different panes.
   int *pane_ids, npanes;
   COM_get_panes( wname.c_str(), &npanes, &pane_ids);
   
   for ( int i=0; i<npanes; ++i) {
     double *coor;
     double *spds;
     COM_get_array( (wname+".nc").c_str(), pane_ids[i], &coor);
     COM_get_array( (wname+".spds").c_str(), pane_ids[i], &spds);
 
     int n;
     COM_get_size( (wname+".nc").c_str(), pane_ids[i], &n);
     init_speed_arrays( coor, spds, n);
   }
 
   COM_free_buffer( &pane_ids);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void init_speed_arrays	(	const double *	coors,
		double *	velocities,
		int	n
	)

Definition at line 361 of file rotate.C.

References i, n, pi, sin, and square().

Referenced by init_speed().

                                                                         {
   const double pi = 3.14159265358979;
   for ( int i=0; i<n; ++i) {
 #if 1
     velocities[3*i] = pi*(-coors[3*i+1]);
     velocities[3*i+1] = pi*(coors[3*i]-1.5);
     velocities[3*i+2] = 0;
 #else
     double xs[]={pi*coors[3*i], pi*coors[3*i+1], pi*coors[3*i+2]};
     velocities[3*i] = 2*square(sin(xs[0]))* sin(2*xs[1])*sin(2*xs[2]);
     velocities[3*i+1] = -sin(2*xs[0])* square(sin(xs[1]))* sin(2*xs[2]);
     velocities[3*i+2] = -sin(2*xs[0])* sin(2*xs[1])* square(sin(xs[2]));
 #endif
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void load_modules ( )

Definition at line 49 of file rotate.C.

References COM_LOAD_MODULE_STATIC_DYNAMIC.

                     {
   COM_LOAD_MODULE_STATIC_DYNAMIC(Rocblas, "BLAS");
   COM_LOAD_MODULE_STATIC_DYNAMIC(Rocmap, "MAP");
   COM_LOAD_MODULE_STATIC_DYNAMIC(Rocprop, "PROP");
   COM_LOAD_MODULE_STATIC_DYNAMIC(Rocsurf, "SURF");
   COM_LOAD_MODULE_STATIC_DYNAMIC(Rocout, "OUT");
 }

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 416 of file rotate.C.

References COM_call_function(), COM_finalize(), COM_get_attribute_handle(), COM_get_function_handle(), COM_init(), COM_print_profile(), COM_set_profiling(), COM_set_profiling_barrier(), COMMPI_Comm_rank(), COMMPI_Initialized(), i, init_attributes(), init_constraints_acmflu(), init_constraints_acmfrac(), init_parameters(), init_speed(), Control_parameter::interval, load_modules(), MPI_COMM_WORLD, output_solution(), Control_parameter::perturb, print_usage(), rank, read_control_file(), read_in_mesh(), Control_parameter::steps, and Control_parameter::timestep.

                                  {
   COM_init( &argc, &argv);
   load_modules();
   print_usage( argc, argv);
 
   if ( COMMPI_Initialized()) rank = COMMPI_Comm_rank( MPI_COMM_WORLD);
 
   // Read in mesh file.
   string wname = read_in_mesh( argv[1]);
 
   // Read in control parameters
   Control_parameter cntr_param;
   read_control_file( argc>2?argv[2]:"control.txt", cntr_param);
   
   // Initialize the attributes and output the initial solution.
   init_attributes( wname, cntr_param);
 
   // Attributes
   int pmesh = COM_get_attribute_handle( (wname+".pmesh").c_str());
   int nc = COM_get_attribute_handle( (wname+".nc").c_str());
   int spds = COM_get_attribute_handle( (wname+".spds").c_str());
   int disps = COM_get_attribute_handle( (wname+".disps").c_str());
   int disps_total = COM_get_attribute_handle( (wname+".disps_total").c_str());
 
   // Funcitions
   int PROP_init = COM_get_function_handle( "PROP.initialize");
   int PROP_perturb = COM_get_function_handle( "PROP.perturb_mesh");
   int PROP_propagate = COM_get_function_handle( "PROP.propagate");
   int BLAS_add = COM_get_function_handle( "BLAS.add");
 
   if ( rank==0) cout << "Propagating interface..." << endl;
   // Initialize control parameters
   COM_call_function( PROP_init, &pmesh);
   init_parameters( cntr_param);
 
   if ( wname.substr(0,6) == "acmflu") 
     init_constraints_acmflu( wname);
   else if ( wname.substr(0,3) == "acm") 
     init_constraints_acmfrac( wname);
   else if ( wname.substr(0,4) == "star") 
     ; // init_constraints_starslice( wname);
 
   if ( cntr_param.perturb > 0)
     COM_call_function( PROP_perturb, &pmesh, &cntr_param.perturb);
 
   output_solution( wname, "00000");
 
   COM_set_profiling(1);
   COM_set_profiling_barrier( PROP_propagate, MPI_COMM_WORLD);
 
   char fname[40];
   std::sprintf(fname, "timedata_%03d.txt", rank);
 
   for ( int i=1; i<=cntr_param.steps; ++i) {
     if ( rank==0) cout << "Step " << i << endl;
     
     init_speed( wname);
     COM_call_function( PROP_propagate, &pmesh, &spds, 
                        &cntr_param.timestep, &disps);
 
     COM_call_function( BLAS_add, &nc, &disps, &nc);
     COM_call_function( BLAS_add, &disps_total, &disps, &disps_total);
 
     if ( i%cntr_param.interval == 0) {
       char   steps[10];
       std::sprintf( steps, "%05d", i);
       output_solution( wname, steps);
 
       COM_print_profile( fname, "Proptest");
     }
   }
 
   COM_finalize();
 }

Here is the call graph for this function:

void output_solution	(	const string &	wname,
		const char *	timelevel
	)

Definition at line 396 of file rotate.C.

References COM_call_function(), COM_get_attribute_handle(), COM_get_function_handle(), and COMMPI_Initialized().

                                                                   {
   static int OUT_write = 0, hdl;
 
   if ( OUT_write==0) {
     OUT_write = COM_get_function_handle( "OUT.write_attribute");
     hdl = COM_get_attribute_handle( (wname+".all").c_str());
 
     int OUT_set = COM_get_function_handle( "OUT.set_option");
     COM_call_function( OUT_set, "format", "HDF");
   }
 
   std::string fname = wname+"_"+timelevel;
 
   if ( !COMMPI_Initialized()) fname.append(".hdf");
   else fname.append("_");
 
   COM_call_function( OUT_write, (char*)fname.c_str(),
                      &hdl, (char*)wname.c_str(), timelevel);
 }

Here is the call graph for this function:

void print_usage	(	int	argc,
		char *	argv[]
	)

Definition at line 59 of file rotate.C.

                                           {
   if ( argc <= 2) {
     cout << "Usage: " << argv[0] << " <surffile> <controlfile>" << std::endl;
     
     exit(-1);
   }
 }

void read_control_file	(	const char *	fname,
		Control_parameter &	cp
	)

Definition at line 89 of file rotate.C.

References COM_assertion_msg, Control_parameter::courant, Control_parameter::eigthres, Control_parameter::fangle, Control_parameter::interval, Control_parameter::method, Control_parameter::normaldif, Control_parameter::perturb, rank, Control_parameter::rediter, Control_parameter::smoother, Control_parameter::speed, Control_parameter::sploc, Control_parameter::steps, Control_parameter::timestep, Control_parameter::verbose, and Control_parameter::wavefrontal.

                                                                   {
   /* Sample control file:
    * method: fo          # method: fo and mp
    * wavefrontal: 1        # wavefrontal condition
    * normaldif: 1        # normal diffusion
    * eigthres: 1.e-4     # threshold for null space: 0..1 (1.e-4 default)
    * courant: 0.5        # courant constant
    * fangle: 35          # feature edge angle: between 0 and 180
    * smoother: angle     # type of mesh-smoothing algorithm
    * rediter: 1          # Number of iterations for vertex redistribution
    * speed: 0.1          # Speed
    * sploc: e            # location of speed: n or e
    * timestep: 0.001     # time step
    * steps: 100          # number of time steps
    * interval: 10        # output intervals
    * verbose: 1          # verbose level
    */
   ifstream is(fname); COM_assertion_msg( is, "File does not exist");
 
   while ( !is.eof()) {
     char buf[255];
     is.get( buf, 255, ':');
     if ( buf[0] == '\0') { is.getline( buf, 255); continue; }
 
     istringstream istr(buf); 
     string keyword; istr >> keyword;
     is.getline( buf, 255, ':'); 
 
     if ( keyword == "method")
       is >> cp.method;
     else if ( keyword == "wavefrontal")
       is >> cp.wavefrontal;
     else if ( keyword == "normaldif")
       is >> cp.normaldif;
     else if ( keyword == "eigthres")
       is >> cp.eigthres;
     else if ( keyword == "courant")
       is >> cp.courant;
     else if ( keyword == "fangle")
       is >> cp.fangle;
     else if ( keyword == "smoother")
       is >> cp.smoother;
     else if ( keyword == "rediter")
       is >> cp.rediter;
     else if ( keyword == "perturb")
       is >> cp.perturb;
     else if ( keyword == "speed")
       is >> cp.speed;
     else if ( keyword == "sploc")
       is >> cp.sploc;
     else if ( keyword == "timestep")
       is >> cp.timestep;
     else if ( keyword == "steps")
       is >> cp.steps;
     else if ( keyword == "interval")
       is >> cp.interval;
     else if ( keyword == "verbose")
       is >> cp.verbose;
     else
       std::cerr << "Unknow keyword " << keyword << std::endl;
     is.getline( buf, 255); 
   }
 
   if ( rank==0) std::cout << " speed is " << cp.speed << std::endl;
 }

std::string read_in_mesh ( const char * fname )

Definition at line 208 of file rotate.C.

References COM_assertion_msg, rank, and IM_Reader::read_winmesh().

                                             {
   if ( rank==0) cout << "Reading surface mesh file \"" << fname << '"' << endl;
 
   std::string fname_str(fname);
 
   std::string::size_type pos = fname_str.find_first_of( ".");
   const string wname = fname_str.substr( 0, pos);
 
   if ( rank==0) cout << "Creating window \"" << wname << '"' << endl;
 
   IM_Reader im_reader;
   int npanes = im_reader.read_winmesh( fname, wname); 
   COM_assertion_msg( npanes>=0, "Failed to read in mesh file. File empty?");
 
   return wname;
 }