Index: source/computed_field/computed_field_matrix_operations.cpp =================================================================== --- source/computed_field/computed_field_matrix_operations.cpp (revision 5305) +++ source/computed_field/computed_field_matrix_operations.cpp (working copy) @@ -1,7 +1,7 @@ /******************************************************************************* FILE : computed_field_matrix_operations.c -LAST MODIFIED : 25 August 2006 +LAST MODIFIED : 24 July 2007 DESCRIPTION : Implements a number of basic vector operations on computed fields. @@ -3017,6 +3017,673 @@ return (return_code); } /* define_Computed_field_type_transpose */ +namespace { + +char computed_field_general_rotation_type_string[] = "general_rotation"; + +class Computed_field_general_rotation : public Computed_field_core +{ +public: + /* cache the rotation matrix */ + double *m; + + Computed_field_general_rotation(Computed_field *field) : + Computed_field_core(field) + { + m = (double*)NULL; + }; + + ~Computed_field_general_rotation(); + +private: + Computed_field_core *copy(Computed_field* new_parent) + { + return new Computed_field_general_rotation(new_parent); + } + + char *get_type_string() + { + return(computed_field_general_rotation_type_string); + } + + int compare(Computed_field_core* other_field); + + int evaluate_cache_at_location(Field_location* location); + + int list(); + + char* get_command_string(); + + int clear_cache(); + int evaluate(); +}; + +Computed_field_general_rotation::~Computed_field_general_rotation() +/******************************************************************************* +LAST MODIFIED : 23 July 2007 + +DESCRIPTION : +Clear the type specific data used by this type. +==============================================================================*/ +{ + ENTER(Computed_field_general_rotation::~Computed_field_general_rotation); + if (field) + { + if (m) + { + DEALLOCATE(m); + } + } + else + { + display_message(ERROR_MESSAGE, + "Computed_field_general_rotation::~Computed_field_general_rotation. " + "Invalid argument(s)"); + } + LEAVE; +} /* Computed_field_general_rotation::~Computed_field_general_rotation */ + +int Computed_field_general_rotation::clear_cache() +/******************************************************************************* +LAST MODIFIED : 23 July 2007 + +DESCRIPTION : +==============================================================================*/ +{ + int return_code; + ENTER(Computed_field_general_rotation::clear_cache); + if (field) + { + if (m) + { + DEALLOCATE(m); + } + return_code = 1; + } + else + { + display_message(ERROR_MESSAGE, + "Computed_field_general_rotation::clear_cache. " + "Invalid argument(s)"); + return_code = 0; + } + LEAVE; + return (return_code); +} /* Computed_field_general_rotation::clear_cache */ + +int Computed_field_general_rotation::compare(Computed_field_core *other_core) +/******************************************************************************* +LAST MODIFIED : 09 July 2007 + +DESCRIPTION : +Compare the type specific data +==============================================================================*/ +{ + Computed_field_general_rotation* other; + int return_code; + + ENTER(Computed_field_general_rotation::compare); + if (field && + (other = dynamic_cast(other_core))) + { + /* no type specific data? */ + return_code = 1; + } + else + { + return_code = 0; + } + LEAVE; + + return (return_code); +} /* Computed_field_general_rotation::compare */ + +int Computed_field_general_rotation::evaluate() +/******************************************************************************* +LAST MODIFIED : 23 July 2007 + +DESCRIPTION : +Evaluates the rotation matrix. The cache is +used to store the matrix in double +precision. Expects that source fields have +been pre-calculated before calling this. +==============================================================================*/ +{ + double rad,u,v,w,L,u2,v2,w2,L2,ca,sa,p1[3]; + int return_code = 0; + + ENTER(Computed_field_general_rotation::evaluate); + if (field) + { + if (m) + { + /* rotation matrix already exists, do nothing? */ + return_code = 1; + } + else if (ALLOCATE(m,double,16)) + { + /* calculate the rotation matrix */ + p1[0] = (double)(field->source_fields[0]->values[0]); + p1[1] = (double)(field->source_fields[0]->values[1]); + p1[2] = (double)(field->source_fields[0]->values[2]); + u = (double)(field->source_fields[1]->values[0]); + v = (double)(field->source_fields[1]->values[1]); + w = (double)(field->source_fields[1]->values[2]); + u2 = u*u; + v2 = v*v; + w2 = w*w; + L2 = u2+v2+w2; + L = sqrt(L2); + rad = (double)(field->source_fields[2]->values[0])*(M_PI/180.0); + ca = cos(rad); + sa = sin(rad); + m[0] = (u2+(v2+w2)*ca)/L2; + m[1] = (u*v*(1-ca)+w*L*sa)/L2; + m[2] = (u*w*(1-ca)-v*L*sa)/L2; + m[3] = 0; + m[4] = (u*v*(1-ca)-w*L*sa)/L2; + m[5] = (v2+(u2+w2)*ca)/L2; + m[6] = (v*w*(1-ca)+u*L*sa)/L2; + m[7] = 0; + m[8] = (u*w*(1-ca)+v*L*sa)/L2; + m[9] = (v*w*(1-ca)-u*L*sa)/L2; + m[10] = (w2+(u2+v2)*ca)/L2; + m[11] = 0; + m[12] = (p1[0]*(v2+w2)-u*(p1[1]*v+p1[2]*w)+(u*(p1[1]*v+p1[2]*w)- + p1[0]*(v2+w2))*ca+(p1[1]*w-p1[2]*v)*L*sa)/L2; + m[13] = (p1[1]*(u2+w2)-v*(p1[0]*u+p1[2]*w)+(v*(p1[0]*u+p1[2]*w)- + p1[1]*(u2+w2))*ca+(p1[2]*u-p1[0]*w)*L*sa)/L2; + m[14] = (p1[2]*(u2+v2)-w*(p1[0]*u+p1[1]*v)+(w*(p1[0]*u+p1[1]*v)- + p1[2]*(u2+v2))*ca+(p1[0]*v-p1[1]*u)*L*sa)/L2; + m[15] = 1; + return_code = 1; + } + else + { + display_message(ERROR_MESSAGE, + "Computed_field_general_rotation::evaluate. Could not allocate cache"); + return_code=0; + } + } + else + { + display_message(ERROR_MESSAGE, + "Computed_field_general_rotation::evaluate. Invalid argument(s)"); + return_code=0; + } + + return (return_code); +} /* Computed_field_general_rotation::evaluate */ + +int Computed_field_general_rotation::evaluate_cache_at_location( + Field_location* location) +/******************************************************************************* +LAST MODIFIED : 09 July 2007 + +DESCRIPTION : +Evaluate the fields cache in the element. +==============================================================================*/ +{ + int i,j,return_code; + //FE_value *temp, *temp2; + FE_value *source,*destination; + int number_of_derivatives,di; + + ENTER(Computed_field_general_rotation::evaluate_cache_at_location); + if (field && location && (field->number_of_source_fields == 4) && + (field->number_of_components >= 3) && + (field->source_fields[0]->number_of_components >= 3) && + (field->source_fields[1]->number_of_components >= 3) && + (field->source_fields[2]->number_of_components >= 1) && + (field->source_fields[3]->number_of_components >= 3)) + { + /* 1. Precalculate any source fields that this field depends on */ + if (return_code = + Computed_field_evaluate_source_fields_cache_at_location(field,location) + && evaluate()) + { + /* 2. Calculate the field */ + /* initialise derivatives */ + if (field->source_fields[3]->derivatives_valid) + { + field->derivatives_valid = 1; + number_of_derivatives = location->get_number_of_derivatives(); + for (di=0;dinumber_of_components);di++) + { + field->derivatives[di] = (FE_value)0.0; + } + } + else + { + /* only worry about derivatives if the source field has + valid derivatives? */ + field->derivatives_valid = 0; + number_of_derivatives = 0; + } + for (i=0;i<3;i++) + { + field->values[i] = (FE_value)0.0; + for (j=0;j<3;j++) + { + /* the coordinate transformation */ + field->values[i] += (FE_value)(m[i+j*4]) * + field->source_fields[3]->values[j]; + if (field->derivatives_valid) + { + source = field->source_fields[3]->derivatives + + j*number_of_derivatives; + destination = field->derivatives + i*number_of_derivatives; + for (di=0;divalues[i] += (FE_value)(m[i+j*4]); + } + } + } + else + { + display_message(ERROR_MESSAGE, + "Computed_field_general_rotation::evaluate_cache_at_location. " + "Invalid argument(s)"); + return_code = 0; + } + LEAVE; + return (return_code); +} /* Computed_field_general_rotation::evaluate_cache_at_location */ + + +int Computed_field_general_rotation::list() +/******************************************************************************* +LAST MODIFIED : 09 July 2007 + +DESCRIPTION : +==============================================================================*/ +{ + int return_code; + + ENTER(List_Computed_field_general_rotation); + if (field) + { + display_message(INFORMATION_MESSAGE, + " origin field : %s\n",field->source_fields[0]->name); + display_message(INFORMATION_MESSAGE, + " axis field : %s\n",field->source_fields[1]->name); + display_message(INFORMATION_MESSAGE, + " angle field : %s\n",field->source_fields[2]->name); + display_message(INFORMATION_MESSAGE, + " source field : %s\n",field->source_fields[3]->name); + return_code = 1; + } + else + { + display_message(ERROR_MESSAGE, + "list_Computed_field_general_rotation. Invalid argument(s)"); + return_code = 0; + } + LEAVE; + + return (return_code); +} /* list_Computed_field_general_rotation */ + +char *Computed_field_general_rotation::get_command_string() +/******************************************************************************* +LAST MODIFIED : 09 July 2007 + +DESCRIPTION : +Returns allocated command string for reproducing field. Includes type. +==============================================================================*/ +{ + char *command_string, *field_name; + int error; + + ENTER(Computed_field_general_rotation::get_command_string); + command_string = (char *)NULL; + if (field) + { + error = 0; + append_string(&command_string, + computed_field_general_rotation_type_string, &error); + append_string(&command_string, " origin ", &error); + if (GET_NAME(Computed_field)(field->source_fields[0], &field_name)) + { + make_valid_token(&field_name); + append_string(&command_string, field_name, &error); + DEALLOCATE(field_name); + } + append_string(&command_string, " axis ", &error); + if (GET_NAME(Computed_field)(field->source_fields[1], &field_name)) + { + make_valid_token(&field_name); + append_string(&command_string, field_name, &error); + DEALLOCATE(field_name); + } + append_string(&command_string, " angle ", &error); + if (GET_NAME(Computed_field)(field->source_fields[2], &field_name)) + { + make_valid_token(&field_name); + append_string(&command_string, field_name, &error); + DEALLOCATE(field_name); + } + append_string(&command_string, " field ", &error); + if (GET_NAME(Computed_field)(field->source_fields[3], &field_name)) + { + make_valid_token(&field_name); + append_string(&command_string, field_name, &error); + DEALLOCATE(field_name); + } + } + else + { + display_message(ERROR_MESSAGE, + "Computed_field_general_rotation::get_command_string. Invalid field"); + } + LEAVE; + + return (command_string); +} /* Computed_field_general_rotation::get_command_string */ + +} //namespace + +int Computed_field_set_type_general_rotation(struct Computed_field *field, + struct Computed_field *origin_field,struct Computed_field *axis_field, + struct Computed_field *angle_field,struct Computed_field *source_field) +/******************************************************************************* +LAST MODIFIED : 09 July 2007 + +DESCRIPTION : +Converts to type COMPUTED_FIELD_GENERAL_ROTATION producing the result +of being rotated degrees about an arbitrary +units . +If function fails, field is guaranteed to be unchanged from its original state, +although its cache may be lost. +==============================================================================*/ +{ + int number_of_source_fields, return_code; + struct Computed_field **temp_source_fields; + + ENTER(Computed_field_set_type_general_rotation); + if (field && source_field && axis_field && angle_field) + { + if ((origin_field->number_of_components >= 3) && + (source_field->number_of_components >= 3) && + (axis_field->number_of_components >= 3)) + { + /* 1. make dynamic allocations for any new type-specific data */ + number_of_source_fields=4; + if (ALLOCATE(temp_source_fields,struct Computed_field *, + number_of_source_fields)) + { + /* 2. free current type-specific data */ + Computed_field_clear_type(field); + /* 3. establish the new type */ + field->number_of_components = source_field->number_of_components; + temp_source_fields[0]=ACCESS(Computed_field)(origin_field); + temp_source_fields[1]=ACCESS(Computed_field)(axis_field); + temp_source_fields[2]=ACCESS(Computed_field)(angle_field); + temp_source_fields[3]=ACCESS(Computed_field)(source_field); + field->source_fields=temp_source_fields; + field->number_of_source_fields=number_of_source_fields; + field->core = new Computed_field_general_rotation(field); + return_code=1; + } + else + { + display_message(ERROR_MESSAGE, + "Computed_field_set_type_general_rotation. " + "Unable to allocate memory"); + DEALLOCATE(temp_source_fields); + return_code=0; + } + } + else + { + display_message(ERROR_MESSAGE, + "Computed_field_set_type_general_rotation. " + "Only 3-D rotations coded"); + return_code=0; + } + } + else + { + display_message(ERROR_MESSAGE, + "Computed_field_set_type_general_rotation. Invalid argument(s)"); + return_code=0; + } + LEAVE; + + return (return_code); +} /* Computed_field_set_type_general_rotation */ + +int Computed_field_get_type_general_rotation(struct Computed_field *field, + struct Computed_field **origin_field,struct Computed_field **axis_field, + struct Computed_field **angle_field,struct Computed_field **source_field) +/******************************************************************************* +LAST MODIFIED : 09 July 2007 + +DESCRIPTION : +If the field is of type COMPUTED_FIELD_GENERAL_ROTATION, the + used by it are returned. +==============================================================================*/ +{ + Computed_field_general_rotation* core; + int return_code; + + ENTER(Computed_field_get_type_general_rotation); + if (field && + (core=dynamic_cast(field->core)) && + origin_field && axis_field && angle_field && source_field) + { + *origin_field = field->source_fields[0]; + *axis_field = field->source_fields[1]; + *angle_field = field->source_fields[2]; + *source_field = field->source_fields[3]; + return_code=1; + } + else + { + display_message(ERROR_MESSAGE, + "Computed_field_get_type_general_rotation. Invalid argument(s)"); + return_code=0; + } + LEAVE; + + return (return_code); +} /* Computed_field_get_type_general_rotation */ + +int define_Computed_field_type_general_rotation(struct Parse_state *state, + void *field_void,void *computed_field_matrix_operations_package_void) +/******************************************************************************* +LAST MODIFIED : 09 July 2007 + +DESCRIPTION : +Converts into type COMPUTED_FIELD_GENERAL_ROTATION (if it is not +already) and allows its contents to be modified. +==============================================================================*/ +{ + char *current_token; + int one=1,return_code; + struct Computed_field *field,*origin,*axis,*angle,*source; + Computed_field_matrix_operations_package + *computed_field_matrix_operations_package; + struct Option_table *option_table; + struct Set_Computed_field_conditional_data set_origin_data,set_axis_data, + set_angle_data,set_source_data; + + ENTER(define_Computed_field_type_general_rotation); + if (state&&(field=(struct Computed_field *)field_void)&& + (computed_field_matrix_operations_package= + (Computed_field_matrix_operations_package *) + computed_field_matrix_operations_package_void)) + { + return_code=1; + /* get valid parameters for general_rotation field */ + origin = (struct Computed_field *)NULL; + axis = (struct Computed_field *)NULL; + angle = (struct Computed_field *)NULL; + source = (struct Computed_field *)NULL; + if (computed_field_general_rotation_type_string == + Computed_field_get_type_string(field)) + { + return_code=Computed_field_get_type_general_rotation(field, + &origin,&axis,&angle,&source); + } + if (return_code) + { + /* ACCESS the source fields for set functions */ + if (origin) + { + ACCESS(Computed_field)(origin); + } + if (axis) + { + ACCESS(Computed_field)(axis); + } + if (angle) + { + ACCESS(Computed_field)(angle); + } + if (source) + { + ACCESS(Computed_field)(source); + } + /* try to handle help first */ + if (current_token=state->current_token) + { + if (!(strcmp(PARSER_HELP_STRING,current_token)&& + strcmp(PARSER_RECURSIVE_HELP_STRING,current_token))) + { + option_table = CREATE(Option_table)(); + Option_table_add_help(option_table, + "A general_rotation field rotates a given coordinate field " + "about an arbitrary axis vector with the given origin. The " + "direction of the vector " + "determines the direction of rotation for positive and " + "negative angles. The angle field is expected to be a single " + "component field with the angle of rotation in degrees."); + /* origin */ + set_origin_data.computed_field_manager= + computed_field_matrix_operations_package->computed_field_manager; + set_origin_data.conditional_function= + Computed_field_has_up_to_3_numerical_components; + set_origin_data.conditional_function_user_data=(void *)NULL; + Option_table_add_entry(option_table,"origin",&origin, + &set_origin_data,set_Computed_field_conditional); + /* axis */ + set_axis_data.computed_field_manager= + computed_field_matrix_operations_package->computed_field_manager; + set_axis_data.conditional_function= + Computed_field_has_up_to_3_numerical_components; + set_axis_data.conditional_function_user_data=(void *)NULL; + Option_table_add_entry(option_table,"axis",&axis, + &set_axis_data,set_Computed_field_conditional); + /* angle */ + set_angle_data.computed_field_manager= + computed_field_matrix_operations_package->computed_field_manager; + set_angle_data.conditional_function= + Computed_field_has_n_components; + set_angle_data.conditional_function_user_data=(void *)(&one); + Option_table_add_entry(option_table,"angle",&angle, + &set_angle_data,set_Computed_field_conditional); + /* source */ + set_source_data.computed_field_manager= + computed_field_matrix_operations_package->computed_field_manager; + set_source_data.conditional_function= + Computed_field_has_numerical_components; + set_source_data.conditional_function_user_data=(void *)(&one); + Option_table_add_entry(option_table,"field",&source, + &set_source_data,set_Computed_field_conditional); + /* parse the option table */ + return_code=Option_table_multi_parse(option_table,state); + DESTROY(Option_table)(&option_table); + } + else + { + option_table = CREATE(Option_table)(); + Option_table_add_help(option_table, + "A general_rotation field rotates a given coordinate field " + "about an arbitrary unit vector. The direction of the vector " + "determines the direction of rotation for positive and " + "negative angles. The angle field is expected to be a single " + "component field with the angle of rotation in degrees."); + /* origin */ + set_origin_data.computed_field_manager= + computed_field_matrix_operations_package->computed_field_manager; + set_origin_data.conditional_function= + Computed_field_has_up_to_3_numerical_components; + set_origin_data.conditional_function_user_data=(void *)NULL; + Option_table_add_entry(option_table,"origin",&origin, + &set_origin_data,set_Computed_field_conditional); + /* axis */ + set_axis_data.computed_field_manager= + computed_field_matrix_operations_package->computed_field_manager; + set_axis_data.conditional_function= + Computed_field_has_up_to_3_numerical_components; + set_axis_data.conditional_function_user_data=(void *)NULL; + Option_table_add_entry(option_table,"axis",&axis, + &set_axis_data,set_Computed_field_conditional); + /* angle */ + set_angle_data.computed_field_manager= + computed_field_matrix_operations_package->computed_field_manager; + set_angle_data.conditional_function= + Computed_field_has_n_components; + set_angle_data.conditional_function_user_data=(void *)(&one); + Option_table_add_entry(option_table,"angle",&angle, + &set_angle_data,set_Computed_field_conditional); + /* source */ + set_source_data.computed_field_manager= + computed_field_matrix_operations_package->computed_field_manager; + set_source_data.conditional_function= + Computed_field_has_numerical_components; + set_source_data.conditional_function_user_data=(void *)(&one); + Option_table_add_entry(option_table,"field",&source, + &set_source_data,set_Computed_field_conditional); + if (return_code=Option_table_multi_parse(option_table,state)) + { + return_code = Computed_field_set_type_general_rotation(field, + origin,axis,angle,source); + } + DESTROY(Option_table)(&option_table); + } + } + else + { + display_message(ERROR_MESSAGE, "Missing command options"); + return_code = 0; + } + if (origin) + { + DEACCESS(Computed_field)(&origin); + } + if (axis) + { + DEACCESS(Computed_field)(&axis); + } + if (angle) + { + DEACCESS(Computed_field)(&angle); + } + if (source) + { + DEACCESS(Computed_field)(&source); + } + } + } + else + { + display_message(ERROR_MESSAGE, + "define_Computed_field_type_general_rotation. Invalid argument(s)"); + return_code=0; + } + LEAVE; + + return (return_code); +} /* define_Computed_field_type_general_rotation */ + int Computed_field_register_types_matrix_operations( struct Computed_field_package *computed_field_package) /******************************************************************************* @@ -3067,6 +3734,11 @@ computed_field_transpose_type_string, define_Computed_field_type_transpose, computed_field_matrix_operations_package); + return_code = + Computed_field_package_add_type(computed_field_package, + computed_field_general_rotation_type_string, + define_Computed_field_type_general_rotation, + computed_field_matrix_operations_package); } else {