54 const int& n_continuously_interpolated_values,
const int& value_id)
58 if ((value_id < -1) || (value_id >= n_continuously_interpolated_values))
60 std::ostringstream error_stream;
61 error_stream <<
"Value_id " << value_id <<
" is out of range."
63 <<
"It can only take the values -1 (position) "
64 <<
"or an integer in the range 0 to "
65 << n_continuously_interpolated_values << std::endl;
68 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
81 const unsigned el_dim =
dim();
84 const unsigned n_shape =
nnode();
93 std::ostringstream error_message;
94 error_message <<
"Dimension mismatch" << std::endl;
95 error_message <<
"The elemental dimension: " <<
dim()
97 <<
" for the jacobian of the mapping to be well-defined"
100 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
105 for (
unsigned i = 0;
i < el_dim;
i++)
108 for (
unsigned j = 0; j < el_dim; j++)
111 jacobian(
i, j) = 0.0;
113 for (
unsigned l = 0; l < n_shape; l++)
115 for (
unsigned k = 0; k < n_shape_type; k++)
137 const unsigned el_dim =
dim();
141 const unsigned n_shape =
nnode();
144 const unsigned n_row =
N2deriv[el_dim];
149 for (
unsigned i = 0;
i < n_row;
i++)
152 for (
unsigned j = 0; j < el_dim; j++)
155 jacobian2(
i, j) = 0.0;
157 for (
unsigned l = 0; l < n_shape; l++)
160 for (
unsigned k = 0; k < n_shape_type; k++)
181 const unsigned n_node =
nnode();
185 const unsigned n_dim_element =
dim();
189 for (
unsigned i = 0;
i < n_dim_element;
i++)
191 for (
unsigned j = 0; j < n_dim_node; j++)
194 interpolated_G(
i, j) = 0.0;
195 for (
unsigned l = 0; l < n_node; l++)
197 for (
unsigned k = 0; k < n_position_type; k++)
199 interpolated_G(
i, j) +=
221 const unsigned el_dim =
dim();
224 const unsigned n_shape =
nnode();
231 std::ostringstream error_message;
232 error_message <<
"Dimension mismatch" << std::endl;
233 error_message <<
"The elemental dimension: " <<
dim()
235 <<
" for the jacobian of the mapping to be well-defined"
238 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
246 for (
unsigned i = 0;
i < el_dim;
i++)
249 for (
unsigned j = 0; j < el_dim; j++)
251 jacobian(
i, j) = 0.0;
252 inverse_jacobian(
i, j) = 0.0;
256 for (
unsigned l = 0; l < n_shape; l++)
259 for (
unsigned k = 0; k < n_shape_type; k++)
269 for (
unsigned i = 0;
i < el_dim;
i++)
271 det *= jacobian(
i,
i);
280 for (
unsigned i = 0;
i < el_dim;
i++)
282 inverse_jacobian(
i,
i) = 1.0 / jacobian(
i,
i);
296 const unsigned n_node =
nnode();
298 for (
unsigned n = 0; n < n_node; n++)
313 const bool& store_local_dof_pt)
316 const unsigned n_node =
nnode();
330 bool hanging_eqn_numbers =
false;
340 std::deque<unsigned long> global_eqn_number_queue;
344 for (
unsigned n = 0; n < n_node; n++)
347 for (
unsigned j = 0; j < n_cont_values; j++)
355 unsigned n_master = hang_info_pt->
nmaster();
357 for (
unsigned m = 0; m < n_master; m++)
364 if (local_eqn_number_done[j][Master_node_pt] ==
false)
369 if (Master_node_pt->
nvalue() < j)
371 std::ostringstream error_stream;
372 error_stream <<
"Master node for " << j
373 <<
"-th value only has "
374 << Master_node_pt->
nvalue() <<
" stored values!"
378 OOMPH_CURRENT_FUNCTION,
379 OOMPH_EXCEPTION_LOCATION);
391 unsigned local_node_index = n_node;
393 for (
unsigned n1 = 0; n1 < n_node; n1++)
397 if (Master_node_pt ==
node_pt(n1))
399 local_node_index = n1;
405 if (local_node_index < n_node)
421 global_eqn_number_queue.push_back(
eqn_number);
423 if (store_local_dof_pt)
441 local_eqn_number_done[j][Master_node_pt] =
true;
443 hanging_eqn_numbers =
true;
454 if (store_local_dof_pt)
461 if (!hanging_eqn_numbers)
473 std::set<Node*> all_nodes;
475 for (
unsigned j = 0; j < n_node; j++)
484 unsigned n_master = hang_info_pt->
nmaster();
487 for (
unsigned m = 0; m < n_master; m++)
491 unsigned old_size = all_nodes.size();
492 all_nodes.insert(master_node_pt);
493 if (all_nodes.size() > old_size)
504 unsigned old_size = all_nodes.size();
505 all_nodes.insert(nod_pt);
506 if (all_nodes.size() > old_size)
524 std::set<std::pair<Data*, unsigned>>& paired_field_data)
528 for (
unsigned n = 0; n < n_internal; n++)
533 const unsigned n_value = dat_pt->
nvalue();
536 for (
unsigned i = 0;
i < n_value;
i++)
538 paired_field_data.insert(std::make_pair(dat_pt,
i));
543 const unsigned n_node = this->
nnode();
544 for (
unsigned n = 0; n < n_node; n++)
549 const unsigned n_value = nod_pt->
nvalue();
552 for (
unsigned i = 0;
i < n_value;
i++)
560 const unsigned nmaster = hang_pt->
nmaster();
563 for (
unsigned m = 0; m < nmaster; m++)
571 paired_field_data.insert(std::make_pair(master_nod_pt,
i));
579 paired_field_data.insert(std::make_pair(nod_pt,
i));
601 unsigned n_nod =
nnode();
604 if (n_nod == 0)
return;
610 unsigned n_dof =
ndof();
621 for (std::map<Node*, unsigned>::iterator it =
627 Node* nod_pt = it->first;
630 unsigned node_number = it->second;
633 for (
unsigned i = 0;
i < dim_nod;
i++)
636 double backup = nod_pt->
x(
i);
640 nod_pt->
x(
i) += eps_fd;
652 for (
unsigned l = 0; l < n_dof; l++)
654 dresidual_dnodal_coordinates(l,
i, node_number) =
655 (res_pls[l] - res[l]) / eps_fd;
660 nod_pt->
x(
i) = backup;
679 const unsigned n_node =
nnode();
694 const unsigned n_dof =
ndof();
703 int local_unknown = 0;
706 for (
unsigned n = 0; n < n_node; n++)
715 for (
unsigned i = 0;
i < n_value;
i++)
723 if (local_unknown >= 0)
726 double*
const value_pt = local_node_pt->
value_pt(
i);
729 const double old_var = *value_pt;
732 *value_pt += fd_step;
743 for (
unsigned m = 0; m < n_dof; m++)
746 jacobian(m, local_unknown) =
747 (newres[m] - residuals[m]) / fd_step;
780 const unsigned n_master = hang_info_pt->
nmaster();
781 for (
unsigned m = 0; m < n_master; m++)
789 if (local_unknown >= 0)
792 double*
const value_pt = master_node_pt->
value_pt(
i);
795 const double old_var = *value_pt;
798 *value_pt += fd_step;
809 for (
unsigned m = 0; m < n_dof; m++)
812 jacobian(m, local_unknown) =
813 (newres[m] - residuals[m]) / fd_step;
859 const unsigned el_dim =
dim();
861 const unsigned n_shape =
nnode();
865 for (
unsigned i = 0;
i < el_dim;
i++)
868 for (
unsigned j = 0; j < el_dim; j++)
871 jacobian(
i, j) = 0.0;
873 for (
unsigned l = 0; l < n_shape; l++)
875 for (
unsigned k = 0; k < n_shape_type; k++)
898 const unsigned el_dim =
dim();
900 const unsigned n_shape =
nnode();
903 const unsigned n_row =
N2deriv[el_dim];
908 for (
unsigned i = 0;
i < n_row;
i++)
911 for (
unsigned j = 0; j < el_dim; j++)
914 jacobian2(
i, j) = 0.0;
916 for (
unsigned l = 0; l < n_shape; l++)
919 for (
unsigned k = 0; k < n_shape_type; k++)
944 const unsigned el_dim =
dim();
946 const unsigned n_shape =
nnode();
953 for (
unsigned i = 0;
i < el_dim;
i++)
956 for (
unsigned j = 0; j < el_dim; j++)
958 jacobian(
i, j) = 0.0;
959 inverse_jacobian(
i, j) = 0.0;
963 for (
unsigned l = 0; l < n_shape; l++)
966 for (
unsigned k = 0; k < n_shape_type; k++)
976 for (
unsigned i = 0;
i < el_dim;
i++)
978 det *= jacobian(
i,
i);
987 for (
unsigned i = 0;
i < el_dim;
i++)
989 inverse_jacobian(
i,
i) = 1.0 / jacobian(
i,
i);
1006 const unsigned n_node =
nnode();
1009 std::set<Data*> all_position_data_pt;
1013 for (
unsigned n = 0; n < n_node; n++)
1022 unsigned n_master = hang_info_pt->
nmaster();
1025 for (
unsigned m = 0; m < n_master; m++)
1031 all_position_data_pt.insert(
1032 dynamic_cast<SolidNode*
>(Master_node_pt)->variable_position_pt());
1039 all_position_data_pt.insert(
1046 return all_position_data_pt.size();
1058 const unsigned n_node =
nnode();
1063 std::set<Data*> all_position_data_pt;
1071 for (
unsigned n = 0; n < n_node; n++)
1080 unsigned n_master = hang_info_pt->
nmaster();
1083 for (
unsigned m = 0; m < n_master; m++)
1090 dynamic_cast<SolidNode*
>(Master_node_pt)->variable_position_pt();
1093 n_old = all_position_data_pt.size();
1094 all_position_data_pt.insert(pos_data_pt);
1097 if (all_position_data_pt.size() > n_old)
1099 all_position_data_vector_pt.push_back(pos_data_pt);
1113 n_old = all_position_data_pt.size();
1114 all_position_data_pt.insert(pos_data_pt);
1117 if (all_position_data_pt.size() > n_old)
1119 all_position_data_vector_pt.push_back(pos_data_pt);
1127 return all_position_data_vector_pt[j];
1138 std::set<Data*>& geometric_data_pt)
1141 const unsigned n_node = this->
nnode();
1142 for (
unsigned j = 0; j < n_node; j++)
1151 unsigned n_master = hang_info_pt->
nmaster();
1154 for (
unsigned m = 0; m < n_master; m++)
1160 geometric_data_pt.insert(
1161 dynamic_cast<SolidNode*
>(Master_node_pt)->variable_position_pt());
1168 geometric_data_pt.insert(
1179 const bool& store_local_dof_pt)
1182 const unsigned n_node =
nnode();
1199 std::map<Node*, bool> local_eqn_number_done;
1205 std::deque<unsigned long> global_eqn_number_queue;
1209 for (
unsigned n = 0; n < n_node; n++)
1218 unsigned n_master = hang_info_pt->
nmaster();
1220 for (
unsigned m = 0; m < n_master; m++)
1227 if (local_eqn_number_done[Master_node_pt] ==
false)
1237 unsigned local_node_index = n_node;
1239 for (
unsigned n1 = 0; n1 < n_node; n1++)
1243 if (Master_node_pt ==
node_pt(n1))
1245 local_node_index = n1;
1251 if (local_node_index < n_node)
1254 for (
unsigned j = 0; j < n_position_type; j++)
1257 for (
unsigned k = 0; k < nodal_dim; k++)
1261 Position_local_eqn_at_node(j, k) =
1270 for (
unsigned j = 0; j < n_position_type; j++)
1273 for (
unsigned k = 0; k < nodal_dim; k++)
1278 ->position_eqn_number(j, k);
1283 global_eqn_number_queue.push_back(
eqn_number);
1285 if (store_local_dof_pt)
1288 &(Master_node_pt->
x_gen(j, k)));
1305 local_eqn_number_done[Master_node_pt] =
true;
1308 Position_local_eqn_at_node;
1319 if (store_local_dof_pt)
1339 const unsigned n_node =
nnode();
1358 const unsigned n_dof =
ndof();
1367 int local_unknown = 0;
1370 for (
unsigned l = 0; l < n_node; l++)
1379 for (
unsigned k = 0; k < n_position_type; k++)
1382 for (
unsigned i = 0;
i < nodal_dim;
i++)
1386 if (local_unknown >= 0)
1389 double*
const value_pt = &(local_node_pt->
x_gen(k,
i));
1392 const double old_var = *value_pt;
1395 *value_pt += fd_step;
1410 for (
unsigned m = 0; m < n_dof; m++)
1413 jacobian(m, local_unknown) =
1414 (newres[m] - residuals[m]) / fd_step;
1437 *value_pt = old_var;
1454 const unsigned n_master = hang_info_pt->
nmaster();
1455 for (
unsigned m = 0; m < n_master; m++)
1465 for (
unsigned k = 0; k < n_position_type; k++)
1468 for (
unsigned i = 0;
i < nodal_dim;
i++)
1470 local_unknown = Position_local_eqn_at_node(k,
i);
1472 if (local_unknown >= 0)
1475 double*
const value_pt = &(master_node_pt->
x_gen(k,
i));
1479 const double old_var = *value_pt;
1482 *value_pt += fd_step;
1496 for (
unsigned m = 0; m < n_dof; m++)
1499 jacobian(m, local_unknown) =
1500 (newres[m] - residuals[m]) / fd_step;
1524 *value_pt = old_var;
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
A class that represents a collection of data; each Data object may contain many different individual ...
static long Is_pinned
Static "Magic number" used in place of the equation number to indicate that the value is pinned.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
double * value_pt(const unsigned &i) const
Return the pointer to the i-the stored value. Typically this is required when direct access to the st...
virtual void update_before_nodal_fd()
Function that is called before the finite differencing of any nodal data. This may be overloaded to u...
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
virtual void update_in_nodal_fd(const unsigned &i)
Function called within the finite difference loop for nodal data after a change in the i-th nodal val...
void check_jacobian(const double &jacobian) const
Helper function used to check for singular or negative Jacobians in the transform from local to globa...
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
double nodal_position_gen(const unsigned &n, const unsigned &k, const unsigned &i) const
Return the value of the k-th type of the i-th positional variable at the local node n.
unsigned nnode() const
Return the number of nodes.
virtual void reset_in_nodal_fd(const unsigned &i)
Function called within the finite difference loop for nodal data after the i-th nodal values is reset...
virtual void reset_after_nodal_fd()
Function that is call after the finite differencing of the nodal data. This may be overloaded to rese...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
static const unsigned N2deriv[]
Static array that holds the number of second derivatives as a function of the dimension of the elemen...
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
unsigned ndof() const
Return the number of equations/dofs in the element.
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
int local_eqn_number(const unsigned long &ieqn_global) const
Return the local equation number corresponding to the ieqn_global-th global equation number....
virtual void get_residuals(Vector< double > &residuals)
Calculate the vector of residuals of the equations in the element. By default initialise the vector t...
unsigned ninternal_data() const
Return the number of internal data objects.
void add_global_eqn_numbers(std::deque< unsigned long > const &global_eqn_numbers, std::deque< double * > const &global_dof_pt)
Add the contents of the queue global_eqn_numbers to the local storage for the local-to-global transla...
static std::deque< double * > Dof_pt_deque
Static storage for deque used to add_global_equation_numbers when pointers to the dofs in each elemen...
Class that contains data for hanging nodes.
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
unsigned nmaster() const
Return the number of master nodes.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
void perform_auxiliary_node_update_fct()
Execute auxiliary update function (if any) – this can be used to update any nodal values following th...
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
bool is_hanging() const
Test whether the node is geometrically hanging.
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i). ‘Type’: k; Coordinate direction: i.
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
void assemble_local_to_eulerian_jacobian2(const DShape &d2psids, DenseMatrix< double > &jacobian2) const
Assemble the the "jacobian" matrix of second derivatives of the mapping from local to Eulerian coordi...
static double Max_integrity_tolerance
Max. allowed discrepancy in element integrity check.
void assign_hanging_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers for hanging node variables.
virtual void deactivate_element()
Final operations that must be performed when the element is no longer active in the mesh,...
double local_to_eulerian_mapping_diagonal(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Eulerian coordinates given the derivatives of the shape functions...
virtual unsigned ncont_interpolated_values() const =0
Number of continuously interpolated values. Note: We assume that they are located at the beginning of...
static void check_value_id(const int &n_continuously_interpolated_values, const int &value_id)
Static helper function that is used to check that the value_id is in range.
void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
virtual void fill_in_jacobian_from_nodal_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the jacobian from the nodal degrees of freedom using finite difference...
void assemble_local_to_eulerian_jacobian(const DShape &dpsids, DenseMatrix< double > &jacobian) const
Assemble the jacobian matrix for the mapping from local to Eulerian coordinates, given the derivative...
std::map< Node *, int > * Local_hang_eqn
Storage for local equation numbers of hanging node variables (values stored at master nodes)....
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Access function that returns the local equation number for the hanging node variables (values stored ...
void assemble_eulerian_base_vectors(const DShape &dpsids, DenseMatrix< double > &interpolated_G) const
Assemble the covariant Eulerian base vectors, assuming that the derivatives of the shape functions wi...
void identify_field_data_for_interactions(std::set< std::pair< Data *, unsigned > > &paired_field_data)
The purpose of this function is to identify all possible Data that can affect the fields interpolated...
std::map< Node *, unsigned > Shape_controlling_node_lookup
Lookup scheme for unique number associated with any of the nodes that actively control the shape of t...
virtual ~RefineableElement()
Destructor, delete the allocated storage for the hanging equations.
Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th Data item that the object's shape depends on: Positional data of non-hangi...
double local_to_lagrangian_mapping_diagonal(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Lagrangian coordinates given the derivatives of the shape functio...
void assemble_local_to_lagrangian_jacobian2(const DShape &d2psids, DenseMatrix< double > &jacobian2) const
Assemble the the "jacobian" matrix of second derivatives, given the second derivatives of the shape f...
unsigned ngeom_data() const
The number of geometric data affecting a RefineableSolidFiniteElement is the positional Data of all n...
void assign_solid_hanging_local_eqn_numbers(const bool &store_local_dof_pt)
Assign local equation numbers to the hanging values associated with positions or additional solid val...
void fill_in_jacobian_from_solid_position_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute element residual Vector and element Jacobian matrix corresponding to the solid positions....
std::map< Node *, DenseMatrix< int > > Local_position_hang_eqn
Storage for local equation numbers of hanging node variables associated with nodal positions....
void identify_geometric_data(std::set< Data * > &geometric_data_pt)
Specify Data that affects the geometry of the element by adding the position Data to the set that's p...
void assemble_local_to_lagrangian_jacobian(const DShape &dpsids, DenseMatrix< double > &jacobian) const
Assemble the jacobian matrix for the mapping from local to lagrangian coordinates,...
double lagrangian_position_gen(const unsigned &n, const unsigned &k, const unsigned &i) const
Return Generalised Lagrangian coordinate at local node n. ‘Direction’ i, ‘Type’ k.
unsigned nnodal_lagrangian_type() const
Return the number of types of (generalised) nodal Lagrangian coordinates required to interpolate the ...
virtual void update_before_solid_position_fd()
Function that is called before the finite differencing of any solid position data....
virtual void update_in_solid_position_fd(const unsigned &i)
Function called within the finite difference loop for the solid position dat after a change in any va...
int position_local_eqn(const unsigned &n, const unsigned &k, const unsigned &j) const
Access function that returns the local equation number that corresponds to the j-th coordinate of the...
virtual void reset_in_solid_position_fd(const unsigned &i)
Function called within the finite difference loop for solid position data after the values in the i-t...
virtual void reset_after_solid_position_fd()
Function that is call after the finite differencing of the solid position data. This may be overloade...
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
A slight extension to the standard template vector class so that we can include "graceful" array rang...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...