78 for (
unsigned i = 0;
i < n_internal_data;
i++)
88 <<
"---------------------------------------------------------------"
91 <<
"Info: Data is already included in this element's internal "
94 <<
"It's stored as entry " <<
i <<
" and I'm not adding it again."
97 <<
"Note: You can suppress this message by recompiling without"
98 <<
"\n PARANOID or setting the boolean \n"
100 "GeneralisedElement::Suppress_warning_about_repeated_internal_"
102 <<
"\n\n to true." << std::endl
103 <<
"---------------------------------------------------------------"
115 Data** new_data_pt =
new Data*[n_internal_data + n_external_data + 1];
118 for (
unsigned i = 0;
i < n_internal_data;
i++)
124 new_data_pt[n_internal_data] = data_pt;
127 for (
unsigned i = 0;
i < n_external_data;
i++)
129 new_data_pt[n_internal_data + 1 +
i] =
Data_pt[n_internal_data +
i];
139 Data_fd.resize(n_internal_data + n_external_data + 1);
141 for (
unsigned i = n_external_data;
i > 0;
i--)
152 return n_internal_data;
162 std::deque<unsigned long>
const& global_eqn_numbers,
163 std::deque<double*>
const& global_dof_pt)
166 const unsigned n_dof =
Ndof;
168 const unsigned n_additional_dof = global_eqn_numbers.size();
170 if (n_additional_dof == 0)
176 const unsigned new_n_dof = n_dof + n_additional_dof;
178 unsigned long* new_eqn_number =
new unsigned long[new_n_dof];
181 for (
unsigned i = 0;
i < n_dof;
i++)
187 unsigned index = n_dof;
189 for (std::deque<unsigned long>::const_iterator it =
190 global_eqn_numbers.begin();
191 it != global_eqn_numbers.end();
195 new_eqn_number[index] = *it;
202 const unsigned n_additional_dof_pt = global_dof_pt.size();
203 if (n_additional_dof_pt > 0)
207 if (n_additional_dof_pt != n_additional_dof)
209 std::ostringstream error_stream;
211 <<
"global_dof_pt is non-empty, yet it does not have the same size\n"
212 <<
"as global_eqn_numbers.\n"
213 <<
"There are " << n_additional_dof <<
" equation numbers,\n"
214 <<
"but " << n_additional_dof_pt << std::endl;
217 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
222 double** new_dof_pt =
new double*[new_n_dof];
224 for (
unsigned i = 0;
i < n_dof;
i++)
230 unsigned index = n_dof;
232 for (std::deque<double*>::const_iterator it = global_dof_pt.begin();
233 it != global_dof_pt.end();
237 new_dof_pt[index] = *it;
281 GeneralisedElement::~GeneralisedElement()
322 for (
unsigned i = 0;
i < n_external_data;
i++)
332 <<
"---------------------------------------------------------------"
335 <<
"Info: Data is already included in this element's external "
338 <<
"It's stored as entry " <<
i <<
" and I'm not adding it again"
341 <<
"Note: You can suppress this message by recompiling without"
342 <<
"\n PARANOID or setting the boolean \n"
344 "GeneralisedElement::Suppress_warning_about_repeated_external_"
346 <<
"\n\n to true." << std::endl
347 <<
"---------------------------------------------------------------"
359 Data** new_data_pt =
new Data*[n_internal_data + n_external_data + 1];
362 for (
unsigned i = 0;
i < (n_internal_data + n_external_data);
i++)
368 new_data_pt[n_internal_data + n_external_data] = data_pt;
377 Data_fd.resize(n_internal_data + n_external_data + 1);
379 Data_fd[n_internal_data + n_external_data] = fd;
385 return n_external_data;
397 if (n_external_data > 0)
400 Data** new_data_pt = 0;
406 if (n_internal_data > 0)
409 new_data_pt =
new Data*[n_internal_data];
411 for (
unsigned i = 0;
i < n_internal_data;
i++)
425 Data_fd.resize(n_internal_data);
441 unsigned index = n_external_data;
443 for (
unsigned i = 0;
i < n_external_data;
i++)
455 if (index < n_external_data)
458 Data** new_data_pt = 0;
464 const unsigned n_total_data = n_internal_data + n_external_data - 1;
467 if (n_total_data > 0)
469 new_data_pt =
new Data*[n_total_data];
473 for (
unsigned i = 0;
i < n_internal_data;
i++)
479 unsigned counter = 0;
480 for (
unsigned i = 0;
i < n_external_data;
i++)
487 new_data_pt[n_internal_data + counter] =
Data_pt[
i];
508 std::ostringstream warning_stream;
510 <<
"Data removed from element's external data " << std::endl
511 <<
"You may have to update the indices for remaining data "
513 <<
"This can be achieved by using add_external_data() "
516 "GeneralisedElement::flush_external_data()",
517 OOMPH_EXCEPTION_LOCATION);
557 std::ostream& out,
const std::string& current_string)
const
561 std::stringstream conversion;
562 conversion <<
" of Internal Data " <<
i << current_string;
579 std::ostream& out,
const std::string& current_string)
const
586 for (
unsigned i = 0;
i < n_internal_data;
i++)
591 std::stringstream conversion;
592 conversion <<
" of Internal Data " <<
i << current_string;
599 for (
unsigned i = 0;
i < n_external_data;
i++)
604 std::stringstream conversion;
605 conversion <<
" of External Data " <<
i << current_string;
617 std::map<unsigned, double*>& map_of_value_pt)
676 const Vector<long>& vector_of_eqn_numbers,
unsigned& index)
697 const bool& store_local_dof_pt)
706 std::ostringstream error_stream;
710 std::map<unsigned, bool> is_repeated;
711 std::set<unsigned long> set_of_global_eqn_numbers;
712 unsigned old_ndof = 0;
713 for (
unsigned n = 0; n <
Ndof; ++n)
715 set_of_global_eqn_numbers.insert(
Eqn_number[n]);
716 if (set_of_global_eqn_numbers.size() == old_ndof)
718 error_stream <<
"Repeated global eqn: " <<
Eqn_number[n] << std::endl;
721 old_ndof = set_of_global_eqn_numbers.size();
725 if (set_of_global_eqn_numbers.size() !=
Ndof)
728 error_stream <<
"Element is ";
729 if (!
is_halo()) error_stream <<
"not ";
730 error_stream <<
"a halo element\n\n";
732 error_stream <<
"\nLocal/lobal equation numbers: " << std::endl;
733 for (
unsigned n = 0; n <
Ndof; ++n)
735 error_stream <<
" " << n <<
" " <<
Eqn_number[n] << std::endl;
739 error_stream << std::endl <<
" element address is " <<
this << std::endl;
743 error_stream <<
"Number of internal data " << nint << std::endl;
744 for (
unsigned i = 0;
i < nint;
i++)
747 unsigned nval = data_pt->
nvalue();
748 for (
unsigned j = 0; j < nval; j++)
751 error_stream <<
"Internal dof: " << eqn_no << std::endl;
752 if (is_repeated[
unsigned(eqn_no)])
754 error_stream <<
"Repeated internal dof: " << eqn_no << std::endl;
761 error_stream <<
"Number of external data " << next << std::endl;
762 for (
unsigned i = 0;
i < next;
i++)
765 unsigned nval = data_pt->
nvalue();
766 for (
unsigned j = 0; j < nval; j++)
769 error_stream <<
"External dof: " << eqn_no << std::endl;
770 if (is_repeated[
unsigned(eqn_no)])
772 error_stream <<
"Repeated external dof: " << eqn_no;
773 Node* nod_pt =
dynamic_cast<Node*
>(data_pt);
776 error_stream <<
" (is a node at: ";
777 unsigned ndim = nod_pt->
ndim();
778 for (
unsigned ii = 0; ii < ndim; ii++)
780 error_stream << nod_pt->
x(ii) <<
" ";
783 error_stream <<
")\n";
798 error_stream <<
"Number of external element data " << next
801 for (
unsigned i = 0;
i < next;
i++)
803 unsigned nval = data_pt[
i]->nvalue();
804 for (
unsigned j = 0; j < nval; j++)
806 int eqn_no = data_pt[
i]->eqn_number(j);
807 error_stream <<
"External element dof: " << eqn_no << std::endl;
808 if (is_repeated[
unsigned(eqn_no)])
810 error_stream <<
"Repeated external element dof: " << eqn_no;
811 Node* nod_pt =
dynamic_cast<Node*
>(data_pt[
i]);
814 error_stream <<
" (is a node at: ";
815 unsigned ndim = nod_pt->
ndim();
816 for (
unsigned ii = 0; ii < ndim; ii++)
818 error_stream << nod_pt->
x(ii) <<
" ";
821 error_stream <<
")\n";
831 error_stream <<
"Number of external element geom data " << next
835 for (
unsigned i = 0;
i < next;
i++)
837 unsigned nval = data_pt[
i]->nvalue();
838 for (
unsigned j = 0; j < nval; j++)
840 int eqn_no = data_pt[
i]->eqn_number(j);
841 error_stream <<
"External element geometric dof: " << eqn_no
843 if (is_repeated[
unsigned(eqn_no)])
846 <<
"Repeated external element geometric dof: " << eqn_no
847 <<
" " << data_pt[
i]->value(j);
848 Node* nod_pt =
dynamic_cast<Node*
>(data_pt[
i]);
851 error_stream <<
" (is a node at: ";
852 unsigned ndim = nod_pt->
ndim();
853 for (
unsigned ii = 0; ii < ndim; ii++)
855 error_stream << nod_pt->
x(
i) <<
" ";
859 error_stream <<
"\n";
870 unsigned n_node = f_el_pt->
nnode();
871 for (
unsigned n = 0; n < n_node; n++)
874 unsigned nval = nod_pt->
nvalue();
875 for (
unsigned j = 0; j < nval; j++)
878 error_stream <<
"Node " << n <<
": Nodal dof: " << eqn_no;
881 if (is_repeated[
unsigned(eqn_no)])
883 error_stream <<
"Node " << n
884 <<
": Repeated nodal dof: " << eqn_no;
887 error_stream <<
" (resized)";
889 error_stream << std::endl;
894 if (solid_nod_pt != 0)
897 unsigned nval = data_pt->
nvalue();
898 for (
unsigned j = 0; j < nval; j++)
901 error_stream <<
"Node " << n <<
": Positional dof: " << eqn_no
903 if (is_repeated[
unsigned(eqn_no)])
905 error_stream <<
"Repeated positional dof: " << eqn_no <<
" "
906 << data_pt->
value(j) << std::endl;
913 n_node = f_el_pt->
nnode();
914 for (
unsigned n = 0; n < n_node; n++)
917 unsigned n_dim = nod_pt->
ndim();
918 error_stream <<
"Node " << n <<
" at ( ";
919 for (
unsigned i = 0;
i < n_dim;
i++)
921 error_stream << nod_pt->
x(
i) <<
" ";
923 error_stream <<
")" << std::endl;
931 error_stream << std::endl << std::endl;
933 <<
"---------------------------------------------------------------"
936 <<
"Note: You can suppress this warning by recompiling without"
937 <<
"\n PARANOID or setting the boolean \n"
939 "GeneralisedElement::Suppress_warning_about_any_repeated_data"
940 <<
"\n\n to true." << std::endl
942 <<
"Only do this if you know what you're doing; repeated equation\n"
943 <<
"numbers are usually a sign of trouble...\n"
944 <<
"---------------------------------------------------------------"
951 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
967 const bool& store_local_dof_pt)
973 const unsigned n_total_data = n_internal_data + n_external_data;
976 if (n_total_data > 0)
986 for (
unsigned i = 1;
i < n_total_data; ++
i)
1000 if (n_total_values == 0)
1015 for (
unsigned i = 0;
i < n_total_values; ++
i)
1021 for (
unsigned i = 1;
i < n_total_data; ++
i)
1032 std::deque<unsigned long> global_eqn_number_queue;
1035 for (
unsigned i = 0;
i < n_internal_data;
i++)
1040 unsigned n_value = data_pt->
nvalue();
1043 for (
unsigned j = 0; j < n_value; j++)
1051 global_eqn_number_queue.push_back(
eqn_number);
1053 if (store_local_dof_pt)
1072 for (
unsigned i = 0;
i < n_external_data;
i++)
1077 unsigned n_value = data_pt->
nvalue();
1080 for (
unsigned j = 0; j < n_value; j++)
1088 global_eqn_number_queue.push_back(
eqn_number);
1090 if (store_local_dof_pt)
1111 if (store_local_dof_pt)
1133 const bool& fd_all_data)
1139 if (n_internal_data == 0)
1149 const unsigned n_dof =
ndof();
1155 int local_unknown = 0;
1161 for (
unsigned i = 0;
i < n_internal_data;
i++)
1171 for (
unsigned j = 0; j < n_value; j++)
1176 if (local_unknown >= 0)
1182 const double old_var = *value_pt;
1185 *value_pt += fd_step;
1194 for (
unsigned m = 0; m < n_dof; m++)
1196 double sum = (newres[m] - residuals[m]) / fd_step;
1198 jacobian(m, local_unknown) = sum;
1202 *value_pt = old_var;
1230 const bool& fd_all_data)
1235 if (n_external_data == 0)
1245 const unsigned n_dof =
ndof();
1251 int local_unknown = 0;
1257 for (
unsigned i = 0;
i < n_external_data;
i++)
1267 for (
unsigned j = 0; j < n_value; j++)
1272 if (local_unknown >= 0)
1278 const double old_var = *value_pt;
1281 *value_pt += fd_step;
1290 for (
unsigned m = 0; m < n_dof; m++)
1292 double sum = (newres[m] - residuals[m]) / fd_step;
1294 jacobian(m, local_unknown) = sum;
1298 *value_pt = old_var;
1324 "Empty fill_in_contribution_to_mass_matrix() has been called.\n";
1326 "This function is called from the default implementation of\n";
1327 error_message +=
"get_mass_matrix();\n";
1328 error_message +=
"and must calculate the residuals vector and mass matrix ";
1329 error_message +=
"without initialising any of their entries.\n\n";
1332 "If you do not wish to use these defaults, you must overload\n";
1333 error_message +=
"get_mass_matrix(), which must initialise the entries\n";
1334 error_message +=
"of the residuals vector and mass matrix to zero.\n";
1338 "GeneralisedElement::fill_in_contribution_to_mass_matrix()",
1339 OOMPH_EXCEPTION_LOCATION);
1356 "Empty fill_in_contribution_to_jacobian_and_mass_matrix() has been ";
1357 error_message +=
"called.\n";
1359 "This function is called from the default implementation of\n";
1360 error_message +=
"get_jacobian_and_mass_matrix();\n";
1362 "and must calculate the residuals vector and mass and jacobian matrices ";
1363 error_message +=
"without initialising any of their entries.\n\n";
1366 "If you do not wish to use these defaults, you must overload\n";
1368 "get_jacobian_and_mass_matrix(), which must initialise the entries\n";
1370 "of the residuals vector, jacobian and mass matrix to zero.\n";
1374 "GeneralisedElement::fill_in_contribution_to_jacobian_and_mass_matrix()",
1375 OOMPH_EXCEPTION_LOCATION);
1390 "Empty fill_in_contribution_to_dresiduals_dparameter() has been ";
1391 error_message +=
"called.\n";
1393 "This function is called from the default implementation of\n";
1394 error_message +=
"get_dresiduals_dparameter();\n";
1395 error_message +=
"and must calculate the derivatives of the residuals "
1396 "vector with respect\n";
1397 error_message +=
"to the passed parameter ";
1398 error_message +=
"without initialising any values.\n\n";
1401 "If you do not wish to use these defaults, you must overload\n";
1403 "get_dresiduals_dparameter(), which must initialise the entries\n";
1404 error_message +=
"of the returned vector to zero.\n";
1407 "This function is intended for use instead of the default (global) \n";
1409 "finite-difference routine when analytic expressions are to be used\n";
1410 error_message +=
"in continuation and bifurcation tracking problems.\n\n";
1411 error_message +=
"This function is only called when the function\n";
1413 "Problem::set_analytic_dparameter() has been called in the driver code\n";
1417 "GeneralisedElement::fill_in_contribution_to_dresiduals_dparameter()",
1418 OOMPH_EXCEPTION_LOCATION);
1431 double*
const& parameter_pt,
1436 "Empty fill_in_contribution_to_djacobian_dparameter() has been ";
1437 error_message +=
"called.\n";
1439 "This function is called from the default implementation of\n";
1440 error_message +=
"get_djacobian_dparameter();\n";
1442 "and must calculate the derivatives of residuals vector and jacobian ";
1443 error_message +=
"matrix\n";
1444 error_message +=
"with respect to the passed parameter ";
1445 error_message +=
"without initialising any values.\n\n";
1448 "If you do not wish to use these defaults, you must overload\n";
1450 "get_djacobian_dparameter(), which must initialise the entries\n";
1451 error_message +=
"of the returned vector and matrix to zero.\n\n";
1454 "This function is intended for use instead of the default (global) \n";
1456 "finite-difference routine when analytic expressions are to be used\n";
1457 error_message +=
"in continuation and bifurcation tracking problems.\n\n";
1458 error_message +=
"This function is only called when the function\n";
1460 "Problem::set_analytic_dparameter() has been called in the driver code\n";
1465 "GeneralisedElement::fill_in_contribution_to_djacobian_dparameter()",
1466 OOMPH_EXCEPTION_LOCATION);
1482 double*
const& parameter_pt,
1488 "Empty fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter() "
1490 error_message +=
" been called.\n";
1492 "This function is called from the default implementation of\n";
1493 error_message +=
"get_djacobian_and_dmass_matrix_dparameter();\n";
1495 "and must calculate the residuals vector and mass and jacobian matrices ";
1496 error_message +=
"without initialising any of their entries.\n\n";
1499 "If you do not wish to use these defaults, you must overload\n";
1500 error_message +=
"get_djacobian_and_dmass_matrix_dparameter(), which must "
1502 error_message +=
"entries of the returned vector and matrices to zero.\n";
1506 "This function is intended for use instead of the default (global) \n";
1508 "finite-difference routine when analytic expressions are to be used\n";
1509 error_message +=
"in continuation and bifurcation tracking problems.\n\n";
1510 error_message +=
"This function is only called when the function\n";
1512 "Problem::set_analytic_dparameter() has been called in the driver code\n";
1516 "GeneralisedElement::fill_in_contribution_to_djacobian_"
1517 "and_dmass_matrix_dparameter()",
1518 OOMPH_EXCEPTION_LOCATION);
1537 "Empty fill_in_contribution_to_hessian_vector_products() has been ";
1538 error_message +=
"called.\n";
1540 "This function is called from the default implementation of\n";
1541 error_message +=
"get_hessian_vector_products(); ";
1542 error_message +=
" and must calculate the products \n";
1543 error_message +=
"of the hessian matrix with the (global) ";
1544 error_message +=
"vectors Y and C\n";
1545 error_message +=
"without initialising any values.\n\n";
1548 "If you do not wish to use these defaults, you must overload\n";
1550 "get_hessian_vector_products(), which must initialise the entries\n";
1551 error_message +=
"of the returned vector to zero.\n\n";
1554 "This function is intended for use instead of the default (global) \n";
1556 "finite-difference routine when analytic expressions are to be used.\n";
1557 error_message +=
"This function is only called when the function\n";
1558 error_message +=
"Problem::set_analytic_hessian_products() has been called "
1559 "in the driver code\n";
1563 "GeneralisedElement::fill_in_contribution_to_hessian_vector_product()",
1564 OOMPH_EXCEPTION_LOCATION);
1572 Vector<std::pair<unsigned, unsigned>>
const& history_index,
1576 "Empty fill_in_contribution_to_inner_products() has been called.\n";
1578 "This function is called from the default implementation of\n";
1579 error_message +=
"get_inner_products();\n ";
1580 error_message +=
"It must calculate the inner products over the element\n";
1581 error_message +=
"of given pairs of history values\n";
1582 error_message +=
"without initialision of the output vector.\n\n";
1585 "If you do not wish to use these defaults, you must overload\n";
1587 "get_inner_products(), which must initialise the entries\n";
1588 error_message +=
"of the returned vector to zero.\n\n";
1591 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1604 "Empty fill_in_contribution_to_inner_product_vectors() has been "
1607 "This function is called from the default implementation of\n";
1608 error_message +=
"get_inner_product_vectors(); ";
1609 error_message +=
" and must calculate the vectors \n";
1610 error_message +=
"corresponding to the input history values\n ";
1612 "that when multiplied by other vectors of history values\n";
1613 error_message +=
"return the appropriate dot products.\n\n";
1614 error_message +=
"The output vectors must not be initialised.\n\n";
1617 "If you do not wish to use these defaults, you must overload\n";
1619 "get_inner_products(), which must initialise the entries\n";
1620 error_message +=
"of the returned vector to zero.\n\n";
1623 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1642 oomph_info <<
"\n ERROR: Failed GeneralisedElement::self_test()!"
1644 oomph_info <<
"for internal data object number: " <<
i << std::endl;
1654 oomph_info <<
"\n ERROR: Failed GeneralisedElement::self_test()!"
1656 oomph_info <<
"for external data object number: " <<
i << std::endl;
1676 namespace Locate_zeta_helpers
1715 const unsigned&
i)
const
1738 std::ostream& out,
const std::string& current_string)
const
1756 std::ostream& out,
const std::string& current_string)
const
1759 const unsigned n_node =
nnode();
1760 for (
unsigned n = 0; n < n_node; n++)
1765 std::stringstream conversion;
1766 conversion <<
" of Node " << n << current_string;
1789 std::ostringstream error_stream;
1791 <<
"Determinant of Jacobian matrix is zero --- "
1792 <<
"singular mapping!\nThe determinant of the "
1793 <<
"jacobian is " << std::fabs(jacobian)
1794 <<
" which is smaller than the treshold "
1796 <<
"You can change this treshold, by specifying "
1797 <<
"FiniteElement::Tolerance_for_singular_jacobian \n"
1798 <<
"Here are the nodal coordinates of the inverted element\n"
1799 <<
"in the form \n\n x,y[,z], hang_status\n\n"
1800 <<
"where hang_status = 1 or 2 for non-hanging or hanging\n"
1801 <<
"nodes respectively (useful for automatic sizing of\n"
1802 <<
"tecplot markers to identify the hanging nodes). \n\n";
1804 unsigned n_nod =
nnode();
1805 unsigned hang_count = 0;
1806 for (
unsigned j = 0; j < n_nod; j++)
1808 for (
unsigned i = 0;
i < n_dim_nod;
i++)
1810 error_stream <<
node_pt(j)->
x(
i) <<
" ";
1814 error_stream <<
" 2";
1819 error_stream <<
" 1";
1821 error_stream << std::endl;
1823 error_stream << std::endl << std::endl;
1827 <<
"NOTE: Offending element is associated with a MacroElement\n"
1828 <<
" AND the element has hanging nodes! \n"
1829 <<
" If an element is thin and highly curved, the \n"
1830 <<
" constraints imposed by\n \n"
1831 <<
" (1) inter-element continuity (imposed by the hanging\n"
1832 <<
" node constraints) and \n\n"
1833 <<
" (2) the need to respect curvilinear domain boundaries\n"
1834 <<
" during mesh refinement (imposed by the element's\n"
1835 <<
" macro element mapping)\n\n"
1836 <<
" may be irreconcilable! \n \n"
1837 <<
" You may have to re-design your base mesh to avoid \n"
1838 <<
" the creation of thin, highly curved elements during\n"
1839 <<
" the refinement process.\n"
1843 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1857 std::ostringstream error_stream;
1858 error_stream <<
"Negative Jacobian in transform from "
1859 <<
"local to global coordinates" << std::endl
1860 <<
" You have an inverted coordinate system!"
1864 <<
"Here are the nodal coordinates of the inverted element\n"
1865 <<
"in the form \n\n x,y[,z], hang_status\n\n"
1866 <<
"where hang_status = 1 or 2 for non-hanging or hanging\n"
1867 <<
"nodes respectively (useful for automatic sizing of\n"
1868 <<
"tecplot markers to identify the hanging nodes). \n\n";
1870 unsigned n_nod =
nnode();
1871 unsigned hang_count = 0;
1872 for (
unsigned j = 0; j < n_nod; j++)
1874 for (
unsigned i = 0;
i < n_dim_nod;
i++)
1876 error_stream <<
node_pt(j)->
x(
i) <<
" ";
1880 error_stream <<
" 2";
1885 error_stream <<
" 1";
1887 error_stream << std::endl;
1889 error_stream << std::endl << std::endl;
1893 <<
"NOTE: The inverted element is associated with a MacroElement\n"
1894 <<
" AND the element has hanging nodes! \n"
1895 <<
" If an element is thin and highly curved, the \n"
1896 <<
" constraints imposed by\n \n"
1897 <<
" (1) inter-element continuity (imposed by the hanging\n"
1898 <<
" node constraints) and \n\n"
1899 <<
" (2) the need to respect curvilinear domain boundaries\n"
1900 <<
" during mesh refinement (imposed by the element's\n"
1901 <<
" macro element mapping)\n\n"
1902 <<
" may be irreconcilable! \n \n"
1903 <<
" You may have to re-design your base mesh to avoid \n"
1904 <<
" the creation of thin, highly curved elements during\n"
1905 <<
" the refinement process.\n"
1912 <<
"If you believe that inverted elements do not cause any\n"
1913 <<
"problems in your specific application you can \n "
1914 <<
"suppress this test by: " << std::endl
1915 <<
" i) setting the (static) flag "
1916 <<
"FiniteElement::Accept_negative_jacobian to be true" << std::endl;
1917 error_stream <<
" ii) switching OFF the PARANOID flag" << std::endl
1923 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1939 const unsigned el_dim =
dim();
1942 const unsigned n_shape =
nnode();
1951 std::ostringstream error_message;
1952 error_message <<
"Dimension mismatch" << std::endl;
1955 <<
" for the jacobian of the mapping to be well-defined"
1958 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1964 for (
unsigned i = 0;
i < el_dim;
i++)
1967 for (
unsigned j = 0; j < el_dim; j++)
1970 jacobian(
i, j) = 0.0;
1972 for (
unsigned l = 0; l < n_shape; l++)
1974 for (
unsigned k = 0; k < n_shape_type; k++)
1997 const unsigned el_dim =
dim();
2001 const unsigned n_shape =
nnode();
2004 const unsigned n_row =
N2deriv[el_dim];
2009 for (
unsigned i = 0;
i < n_row;
i++)
2012 for (
unsigned j = 0; j < el_dim; j++)
2015 jacobian2(
i, j) = 0.0;
2017 for (
unsigned l = 0; l < n_shape; l++)
2020 for (
unsigned k = 0; k < n_shape_type; k++)
2043 const unsigned n_node =
nnode();
2047 const unsigned n_dim_element =
dim();
2051 for (
unsigned i = 0;
i < n_dim_element;
i++)
2053 for (
unsigned j = 0; j < n_dim_node; j++)
2056 interpolated_G(
i, j) = 0.0;
2057 for (
unsigned l = 0; l < n_node; l++)
2059 for (
unsigned k = 0; k < n_position_type; k++)
2061 interpolated_G(
i, j) +=
2074 double FiniteElement::invert_jacobian<0>(
2079 oomph_info <<
"\nWarning: You are trying to invert the jacobian for "
2080 <<
"a 'point' element" << std::endl
2081 <<
"This makes no sense and is almost certainly an error"
2094 double FiniteElement::invert_jacobian<1>(
2099 const double det = jacobian(0, 0);
2107 inverse_jacobian(0, 0) = 1.0 / jacobian(0, 0);
2117 double FiniteElement::invert_jacobian<2>(
2123 jacobian(0, 0) * jacobian(1, 1) - jacobian(0, 1) * jacobian(1, 0);
2131 inverse_jacobian(0, 0) = jacobian(1, 1) / det;
2132 inverse_jacobian(0, 1) = -jacobian(0, 1) / det;
2133 inverse_jacobian(1, 0) = -jacobian(1, 0) / det;
2134 inverse_jacobian(1, 1) = jacobian(0, 0) / det;
2145 double FiniteElement::invert_jacobian<3>(
2150 const double det = jacobian(0, 0) * jacobian(1, 1) * jacobian(2, 2) +
2151 jacobian(0, 1) * jacobian(1, 2) * jacobian(2, 0) +
2152 jacobian(0, 2) * jacobian(1, 0) * jacobian(2, 1) -
2153 jacobian(0, 0) * jacobian(1, 2) * jacobian(2, 1) -
2154 jacobian(0, 1) * jacobian(1, 0) * jacobian(2, 2) -
2155 jacobian(0, 2) * jacobian(1, 1) * jacobian(2, 0);
2163 inverse_jacobian(0, 0) =
2164 (jacobian(1, 1) * jacobian(2, 2) - jacobian(1, 2) * jacobian(2, 1)) / det;
2165 inverse_jacobian(0, 1) =
2166 -(jacobian(0, 1) * jacobian(2, 2) - jacobian(0, 2) * jacobian(2, 1)) /
2168 inverse_jacobian(0, 2) =
2169 (jacobian(0, 1) * jacobian(1, 2) - jacobian(0, 2) * jacobian(1, 1)) / det;
2170 inverse_jacobian(1, 0) =
2171 -(jacobian(1, 0) * jacobian(2, 2) - jacobian(1, 2) * jacobian(2, 0)) /
2173 inverse_jacobian(1, 1) =
2174 (jacobian(0, 0) * jacobian(2, 2) - jacobian(0, 2) * jacobian(2, 0)) / det;
2175 inverse_jacobian(1, 2) =
2176 -(jacobian(0, 0) * jacobian(1, 2) - jacobian(0, 2) * jacobian(1, 0)) /
2178 inverse_jacobian(2, 0) =
2179 (jacobian(1, 0) * jacobian(2, 1) - jacobian(1, 1) * jacobian(2, 0)) / det;
2180 inverse_jacobian(2, 1) =
2181 -(jacobian(0, 0) * jacobian(2, 1) - jacobian(0, 1) * jacobian(2, 0)) /
2183 inverse_jacobian(2, 2) =
2184 (jacobian(0, 0) * jacobian(1, 1) - jacobian(0, 1) * jacobian(1, 0)) / det;
2201 const unsigned el_dim =
dim();
2207 return invert_jacobian<0>(jacobian, inverse_jacobian);
2210 return invert_jacobian<1>(jacobian, inverse_jacobian);
2213 return invert_jacobian<2>(jacobian, inverse_jacobian);
2216 return invert_jacobian<3>(jacobian, inverse_jacobian);
2220 std::ostringstream error_stream;
2221 error_stream <<
"Dimension of the element must be 0,1,2 or 3, not "
2222 << el_dim << std::endl;
2225 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2236 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<0>(
2242 oomph_info <<
"\nWarning: You are trying to calculate derivatives of "
2243 <<
"a jacobian w.r.t. nodal coordinates for a 'point' "
2244 <<
"element." << std::endl
2245 <<
"This makes no sense and is almost certainly an error."
2255 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<1>(
2261 const unsigned n_node =
nnode();
2264 for (
unsigned j = 0; j < n_node; j++)
2266 djacobian_dX(0, j) = dpsids(j, 0);
2275 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<2>(
2281 const unsigned n_node =
nnode();
2284 for (
unsigned j = 0; j < n_node; j++)
2287 djacobian_dX(0, j) =
2288 dpsids(j, 0) * jacobian(1, 1) - dpsids(j, 1) * jacobian(0, 1);
2291 djacobian_dX(1, j) =
2292 dpsids(j, 1) * jacobian(0, 0) - dpsids(j, 0) * jacobian(1, 0);
2301 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<3>(
2307 const unsigned n_node =
nnode();
2310 for (
unsigned j = 0; j < n_node; j++)
2313 djacobian_dX(0, j) =
2315 (jacobian(1, 1) * jacobian(2, 2) - jacobian(1, 2) * jacobian(2, 1)) +
2317 (jacobian(0, 2) * jacobian(2, 1) - jacobian(0, 1) * jacobian(2, 2)) +
2319 (jacobian(0, 1) * jacobian(1, 2) - jacobian(0, 2) * jacobian(1, 1));
2322 djacobian_dX(1, j) =
2324 (jacobian(1, 2) * jacobian(2, 0) - jacobian(1, 0) * jacobian(2, 2)) +
2326 (jacobian(0, 0) * jacobian(2, 2) - jacobian(0, 2) * jacobian(2, 0)) +
2328 (jacobian(0, 2) * jacobian(1, 0) - jacobian(0, 0) * jacobian(1, 2));
2331 djacobian_dX(2, j) =
2333 (jacobian(1, 0) * jacobian(2, 1) - jacobian(1, 1) * jacobian(2, 0)) +
2335 (jacobian(0, 1) * jacobian(2, 0) - jacobian(0, 0) * jacobian(2, 1)) +
2337 (jacobian(0, 0) * jacobian(1, 1) - jacobian(0, 1) * jacobian(1, 0));
2347 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<0>(
2348 const double& det_jacobian,
2356 oomph_info <<
"\nWarning: You are trying to calculate derivatives of "
2357 <<
"eulerian derivatives of shape functions w.r.t. nodal "
2358 <<
"coordinates for a 'point' element." << std::endl
2359 <<
"This makes no sense and is almost certainly an error."
2370 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<1>(
2371 const double& det_jacobian,
2379 const double inv_det_jac = 1.0 / det_jacobian;
2382 const unsigned n_node =
nnode();
2385 for (
unsigned q = 0; q < n_node; q++)
2388 for (
unsigned j = 0; j < n_node; j++)
2390 d_dpsidx_dX(0, q, j, 0) =
2391 -djacobian_dX(0, q) * dpsids(j, 0) * inv_det_jac * inv_det_jac;
2402 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<2>(
2403 const double& det_jacobian,
2411 const double inv_det_jac = 1.0 / det_jacobian;
2414 const unsigned n_node =
nnode();
2417 for (
unsigned p = 0; p < 2; p++)
2420 for (
unsigned q = 0; q < n_node; q++)
2423 for (
unsigned j = 0; j < n_node; j++)
2426 d_dpsidx_dX(p, q, j, 0) =
2427 -djacobian_dX(p, q) * (inverse_jacobian(0, 0) * dpsids(j, 0) +
2428 inverse_jacobian(0, 1) * dpsids(j, 1));
2432 d_dpsidx_dX(p, q, j, 0) +=
2433 dpsids(j, 0) * dpsids(q, 1) - dpsids(j, 1) * dpsids(q, 0);
2435 d_dpsidx_dX(p, q, j, 0) *= inv_det_jac;
2438 d_dpsidx_dX(p, q, j, 1) =
2439 -djacobian_dX(p, q) * (inverse_jacobian(1, 1) * dpsids(j, 1) +
2440 inverse_jacobian(1, 0) * dpsids(j, 0));
2444 d_dpsidx_dX(p, q, j, 1) +=
2445 dpsids(j, 1) * dpsids(q, 0) - dpsids(j, 0) * dpsids(q, 1);
2447 d_dpsidx_dX(p, q, j, 1) *= inv_det_jac;
2459 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<3>(
2460 const double& det_jacobian,
2468 const double inv_det_jac = 1.0 / det_jacobian;
2471 const unsigned n_node =
nnode();
2474 for (
unsigned p = 0; p < 3; p++)
2477 for (
unsigned q = 0; q < n_node; q++)
2480 for (
unsigned j = 0; j < n_node; j++)
2483 for (
unsigned i = 0;
i < 3;
i++)
2485 d_dpsidx_dX(p, q, j,
i) =
2486 -djacobian_dX(p, q) * (inverse_jacobian(
i, 0) * dpsids(j, 0) +
2487 inverse_jacobian(
i, 1) * dpsids(j, 1) +
2488 inverse_jacobian(
i, 2) * dpsids(j, 2));
2495 d_dpsidx_dX(p, q, j, 1) += ((dpsids(q, 2) * jacobian(1, 2) -
2496 dpsids(q, 1) * jacobian(2, 2)) *
2498 (dpsids(q, 0) * jacobian(2, 2) -
2499 dpsids(q, 2) * jacobian(0, 2)) *
2501 (dpsids(q, 1) * jacobian(0, 2) -
2502 dpsids(q, 0) * jacobian(1, 2)) *
2505 d_dpsidx_dX(p, q, j, 2) += ((dpsids(q, 1) * jacobian(2, 1) -
2506 dpsids(q, 2) * jacobian(1, 1)) *
2508 (dpsids(q, 2) * jacobian(0, 1) -
2509 dpsids(q, 0) * jacobian(2, 1)) *
2511 (dpsids(q, 0) * jacobian(1, 1) -
2512 dpsids(q, 1) * jacobian(0, 1)) *
2518 d_dpsidx_dX(p, q, j, 0) += ((dpsids(q, 1) * jacobian(2, 2) -
2519 dpsids(q, 2) * jacobian(1, 2)) *
2521 (dpsids(q, 2) * jacobian(0, 2) -
2522 dpsids(q, 0) * jacobian(2, 2)) *
2524 (dpsids(q, 0) * jacobian(1, 2) -
2525 dpsids(q, 1) * jacobian(0, 2)) *
2528 d_dpsidx_dX(p, q, j, 2) += ((dpsids(q, 2) * jacobian(1, 0) -
2529 dpsids(q, 1) * jacobian(2, 0)) *
2531 (dpsids(q, 0) * jacobian(2, 0) -
2532 dpsids(q, 2) * jacobian(0, 0)) *
2534 (dpsids(q, 1) * jacobian(0, 0) -
2535 dpsids(q, 0) * jacobian(1, 0)) *
2541 d_dpsidx_dX(p, q, j, 0) += ((dpsids(q, 2) * jacobian(1, 1) -
2542 dpsids(q, 1) * jacobian(2, 1)) *
2544 (dpsids(q, 0) * jacobian(2, 1) -
2545 dpsids(q, 2) * jacobian(0, 1)) *
2547 (dpsids(q, 1) * jacobian(0, 1) -
2548 dpsids(q, 0) * jacobian(1, 1)) *
2551 d_dpsidx_dX(p, q, j, 1) += ((dpsids(q, 1) * jacobian(2, 0) -
2552 dpsids(q, 2) * jacobian(1, 0)) *
2554 (dpsids(q, 2) * jacobian(0, 0) -
2555 dpsids(q, 0) * jacobian(2, 0)) *
2557 (dpsids(q, 0) * jacobian(1, 0) -
2558 dpsids(q, 1) * jacobian(0, 0)) *
2564 for (
unsigned i = 0;
i < 3;
i++)
2566 d_dpsidx_dX(p, q, j,
i) *= inv_det_jac;
2624 const unsigned el_dim =
dim();
2627 const unsigned n_shape =
nnode();
2634 std::ostringstream error_message;
2635 error_message <<
"Dimension mismatch" << std::endl;
2638 <<
" for the jacobian of the mapping to be well-defined"
2641 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2649 for (
unsigned i = 0;
i < el_dim;
i++)
2652 for (
unsigned j = 0; j < el_dim; j++)
2654 jacobian(
i, j) = 0.0;
2655 inverse_jacobian(
i, j) = 0.0;
2659 for (
unsigned l = 0; l < n_shape; l++)
2662 for (
unsigned k = 0; k < n_shape_type; k++)
2672 for (
unsigned i = 0;
i < el_dim;
i++)
2674 det *= jacobian(
i,
i);
2683 for (
unsigned i = 0;
i < el_dim;
i++)
2685 inverse_jacobian(
i,
i) = 1.0 / jacobian(
i,
i);
2705 const unsigned el_dim =
dim();
2709 const unsigned n_node =
nnode();
2712 if (djacobian_dX.
nrow() != el_dim)
2714 std::ostringstream error_message;
2715 error_message <<
"djacobian_dX must have the same number of rows as the"
2716 <<
"\nspatial dimension of the element.";
2718 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2721 if (djacobian_dX.
ncol() != n_node)
2723 std::ostringstream error_message;
2725 <<
"djacobian_dX must have the same number of columns as the"
2726 <<
"\nnumber of nodes in the element.";
2728 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2737 dJ_eulerian_dnodal_coordinates_templated_helper<0>(
2738 jacobian, dpsids, djacobian_dX);
2741 dJ_eulerian_dnodal_coordinates_templated_helper<1>(
2742 jacobian, dpsids, djacobian_dX);
2745 dJ_eulerian_dnodal_coordinates_templated_helper<2>(
2746 jacobian, dpsids, djacobian_dX);
2749 dJ_eulerian_dnodal_coordinates_templated_helper<3>(
2750 jacobian, dpsids, djacobian_dX);
2754 std::ostringstream error_stream;
2755 error_stream <<
"Dimension of the element must be 0,1,2 or 3, not "
2756 << el_dim << std::endl;
2759 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2780 const double& det_jacobian,
2788 const unsigned el_dim =
dim();
2792 const unsigned n_node =
nnode();
2795 if (d_dpsidx_dX.
nindex1() != el_dim || d_dpsidx_dX.
nindex2() != n_node ||
2796 d_dpsidx_dX.
nindex3() != n_node || d_dpsidx_dX.
nindex4() != el_dim)
2798 std::ostringstream error_message;
2799 error_message <<
"d_dpsidx_dX must be of the following dimensions:"
2800 <<
"\nd_dpsidx_dX(el_dim,n_node,n_node,el_dim)";
2802 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2811 d_dshape_eulerian_dnodal_coordinates_templated_helper<0>(
2820 d_dshape_eulerian_dnodal_coordinates_templated_helper<1>(
2829 d_dshape_eulerian_dnodal_coordinates_templated_helper<2>(
2838 d_dshape_eulerian_dnodal_coordinates_templated_helper<3>(
2848 std::ostringstream error_stream;
2849 error_stream <<
"Dimension of the element must be 0,1,2 or 3, not "
2850 << el_dim << std::endl;
2853 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2867 const unsigned n_basis = dbasis.
nindex1();
2868 const unsigned n_basis_type = dbasis.
nindex2();
2870 const unsigned n_dim =
dim();
2873 double new_derivatives[n_dim];
2876 for (
unsigned l = 0; l < n_basis; l++)
2879 for (
unsigned k = 0; k < n_basis_type; k++)
2882 for (
unsigned j = 0; j < n_dim; j++)
2885 new_derivatives[j] = 0.0;
2887 for (
unsigned i = 0;
i < n_dim;
i++)
2889 new_derivatives[j] += inverse_jacobian(j,
i) * dbasis(l, k,
i);
2893 for (
unsigned j = 0; j < n_dim; j++)
2895 dbasis(l, k, j) = new_derivatives[j];
2911 const unsigned n_basis = dbasis.
nindex1();
2912 const unsigned n_basis_type = dbasis.
nindex2();
2914 const unsigned n_dim =
dim();
2917 for (
unsigned l = 0; l < n_basis; l++)
2920 for (
unsigned k = 0; k < n_basis_type; k++)
2923 for (
unsigned j = 0; j < n_dim; j++)
2925 dbasis(l, k, j) *= inverse_jacobian(j, j);
2939 void FiniteElement::transform_second_derivatives_template<1>(
2947 const unsigned n_basis = dbasis.
nindex1();
2948 const unsigned n_basis_type = dbasis.
nindex2();
2952 for (
unsigned l = 0; l < n_basis; l++)
2955 for (
unsigned k = 0; k < n_basis_type; k++)
2958 d2basis(l, k, 0) / (jacobian(0, 0) * jacobian(0, 0))
2960 - dbasis(l, k, 0) * jacobian2(0, 0) /
2961 (jacobian(0, 0) * jacobian(0, 0) * jacobian(0, 0));
2978 void FiniteElement::transform_second_derivatives_template<2>(
2986 const unsigned n_basis = dbasis.
nindex1();
2987 const unsigned n_basis_type = dbasis.
nindex2();
2991 jacobian(0, 0) * jacobian(1, 1) - jacobian(0, 1) * jacobian(1, 0);
3000 jacobian2(0, 0) * jacobian(1, 1) + jacobian(0, 0) * jacobian2(2, 1) -
3001 jacobian2(0, 1) * jacobian(1, 0) - jacobian(0, 1) * jacobian2(2, 0);
3003 jacobian2(2, 0) * jacobian(1, 1) + jacobian(0, 0) * jacobian2(1, 1) -
3004 jacobian2(2, 1) * jacobian(1, 0) - jacobian(0, 1) * jacobian2(1, 0);
3008 double dinverse_jacobiands[2][2][2];
3010 dinverse_jacobiands[0][0][0] =
3011 jacobian2(2, 1) / det - jacobian(1, 1) * ddetds[0] / (det * det);
3012 dinverse_jacobiands[0][1][0] =
3013 -jacobian2(0, 1) / det + jacobian(0, 1) * ddetds[0] / (det * det);
3014 dinverse_jacobiands[1][0][0] =
3015 -jacobian2(2, 0) / det + jacobian(1, 0) * ddetds[0] / (det * det);
3016 dinverse_jacobiands[1][1][0] =
3017 jacobian2(0, 0) / det - jacobian(0, 0) * ddetds[0] / (det * det);
3019 dinverse_jacobiands[0][0][1] =
3020 jacobian2(1, 1) / det - jacobian(1, 1) * ddetds[1] / (det * det);
3021 dinverse_jacobiands[0][1][1] =
3022 -jacobian2(2, 1) / det + jacobian(0, 1) * ddetds[1] / (det * det);
3023 dinverse_jacobiands[1][0][1] =
3024 -jacobian2(1, 0) / det + jacobian(1, 0) * ddetds[1] / (det * det);
3025 dinverse_jacobiands[1][1][1] =
3026 jacobian2(2, 0) / det - jacobian(0, 0) * ddetds[1] / (det * det);
3029 DShape dphidXds0(n_basis, n_basis_type, 2),
3030 dphidXds1(n_basis, n_basis_type, 2);
3032 for (
unsigned l = 0; l < n_basis; l++)
3034 for (
unsigned k = 0; k < n_basis_type; k++)
3036 for (
unsigned j = 0; j < 2; j++)
3040 dphidXds0(l, k, j) = dinverse_jacobiands[j][0][0] * dbasis(l, k, 0) +
3041 dinverse_jacobiands[j][1][0] * dbasis(l, k, 1) +
3042 inverse_jacobian(j, 0) * d2basis(l, k, 0) +
3043 inverse_jacobian(j, 1) * d2basis(l, k, 2);
3045 dphidXds1(l, k, j) = dinverse_jacobiands[j][0][1] * dbasis(l, k, 0) +
3046 dinverse_jacobiands[j][1][1] * dbasis(l, k, 1) +
3047 inverse_jacobian(j, 0) * d2basis(l, k, 2) +
3048 inverse_jacobian(j, 1) * d2basis(l, k, 1);
3054 for (
unsigned l = 0; l < n_basis; l++)
3056 for (
unsigned k = 0; k < n_basis_type; k++)
3059 for (
unsigned j = 0; j < 3; j++)
3061 d2basis(l, k, j) = 0.0;
3065 for (
unsigned i = 0;
i < 2;
i++)
3067 d2basis(l, k,
i) = inverse_jacobian(
i, 0) * dphidXds0(l, k,
i) +
3068 inverse_jacobian(
i, 1) * dphidXds1(l, k,
i);
3071 d2basis(l, k, 2) += inverse_jacobian(0, 0) * dphidXds0(l, k, 1) +
3072 inverse_jacobian(0, 1) * dphidXds1(l, k, 1);
3090 void FiniteElement::transform_second_derivatives_diagonal<1>(
3097 FiniteElement::transform_second_derivatives_template<1>(
3098 jacobian, inverse_jacobian, jacobian2, dbasis, d2basis);
3108 void FiniteElement::transform_second_derivatives_diagonal<2>(
3116 const unsigned n_basis = dbasis.
nindex1();
3117 const unsigned n_basis_type = dbasis.
nindex2();
3123 for (
unsigned l = 0; l < n_basis; l++)
3125 for (
unsigned k = 0; k < n_basis_type; k++)
3129 d2basis(l, k, 0) / (jacobian(0, 0) * jacobian(0, 0)) -
3130 dbasis(l, k, 0) * jacobian2(0, 0) /
3131 (jacobian(0, 0) * jacobian(0, 0) * jacobian(0, 0));
3134 d2basis(l, k, 1) / (jacobian(1, 1) * jacobian(1, 1)) -
3135 dbasis(l, k, 1) * jacobian2(1, 1) /
3136 (jacobian(1, 1) * jacobian(1, 1) * jacobian(1, 1));
3138 d2basis(l, k, 2) = d2basis(l, k, 2) / (jacobian(0, 0) * jacobian(1, 1));
3163 const unsigned el_dim =
dim();
3168 transform_second_derivatives_template<1>(
3169 jacobian, inverse_jacobian, jacobian2, dbasis, d2basis);
3172 transform_second_derivatives_template<2>(
3173 jacobian, inverse_jacobian, jacobian2, dbasis, d2basis);
3177 throw OomphLibError(
"Not implemented yet ... maybe one day",
3178 OOMPH_CURRENT_FUNCTION,
3179 OOMPH_EXCEPTION_LOCATION);
3187 std::ostringstream error_stream;
3188 error_stream <<
"Dimension of the element must be 1,2 or 3, not "
3189 << el_dim << std::endl;
3192 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3227 unsigned n_coordinates = s_fraction.size();
3228 for (
unsigned i = 0;
i < n_coordinates;
i++)
3253 const unsigned el_dim =
dim();
3257 for (
unsigned i = 0;
i < el_dim;
i++)
3274 const unsigned el_dim =
dim();
3278 for (
unsigned i = 0;
i < el_dim;
i++)
3310 const unsigned el_dim =
dim();
3314 for (
unsigned i = 0;
i < el_dim;
i++)
3333 const unsigned el_dim =
dim();
3360 const unsigned el_dim =
dim();
3388 const unsigned el_dim =
dim();
3391 unsigned nnod =
nnode();
3392 DShape dpsi(nnod, el_dim);
3421 const unsigned& ipt,
3428 const unsigned el_dim =
dim();
3449 det, jacobian, djacobian_dX, inverse_jacobian, dpsi, d_dpsidx_dX);
3486 const unsigned el_dim =
dim();
3488 const unsigned n_deriv =
N2deriv[el_dim];
3506 jacobian, inverse_jacobian, jacobian2, dpsi, d2psi);
3540 const unsigned el_dim =
dim();
3542 const unsigned n_deriv =
N2deriv[el_dim];
3560 jacobian, inverse_jacobian, jacobian2, dpsi, d2psi);
3575 const bool& store_local_dof_pt)
3578 const unsigned n_node =
nnode();
3589 for (
unsigned n = 1; n < n_node; n++)
3603 if (n_total_values == 0)
3615 for (
unsigned i = 0;
i < n_total_values;
i++)
3621 for (
unsigned n = 1; n < n_node; ++n)
3633 std::deque<unsigned long> global_eqn_number_queue;
3636 for (
unsigned n = 0; n < n_node; n++)
3642 unsigned n_value = nod_pt->
nvalue();
3645 for (
unsigned j = 0; j < n_value; j++)
3653 global_eqn_number_queue.push_back(
eqn_number);
3655 if (store_local_dof_pt)
3676 if (store_local_dof_pt)
3694 const unsigned n_node =
nnode();
3706 const unsigned n_dof =
ndof();
3711 int local_unknown = 0;
3717 for (
unsigned n = 0; n < n_node; n++)
3723 for (
unsigned i = 0;
i < n_value;
i++)
3728 if (local_unknown >= 0)
3734 const double old_var = *value_pt;
3737 *value_pt += fd_step;
3746 for (
unsigned m = 0; m < n_dof; m++)
3748 double sum = (newres[m] - residuals[m]) / fd_step;
3750 jacobian(m, local_unknown) = sum;
3754 *value_pt = old_var;
3778 unsigned n_nod =
nnode();
3781 if (n_nod == 0)
return;
3787 unsigned n_dof =
ndof();
3798 for (
unsigned j = 0; j < n_nod; j++)
3804 for (
unsigned i = 0;
i < dim_nod;
i++)
3807 double backup = nod_pt->
x(
i);
3811 nod_pt->
x(
i) += eps_fd;
3823 for (
unsigned l = 0; l < n_dof; l++)
3825 dresidual_dnodal_coordinates(l,
i, j) =
3826 (res_pls[l] - res[l]) / eps_fd;
3831 nod_pt->
x(
i) = backup;
3849 unsigned n_node =
nnode();
3855 std::vector<int> local_node_number;
3857 for (
unsigned i = 0;
i < n_node;
i++)
3864 local_node_number.push_back(
i);
3871 std::ostringstream error_stream;
3872 error_stream <<
"Node " << global_node_pt <<
" appears " << count
3873 <<
" times in an element." << std::endl
3874 <<
"In positions: ";
3875 for (std::vector<int>::iterator it = local_node_number.begin();
3876 it != local_node_number.end();
3879 error_stream << *it <<
" ";
3881 error_stream << std::endl <<
"That seems very odd." << std::endl;
3884 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3890 for (
unsigned i = 0;
i < n_node;
i++)
3920 const unsigned n_node =
Nnode;
3922 for (
unsigned n = 0; n < n_node; n++)
3927 for (
unsigned i = 0;
i < el_dim;
i++)
3932 if (std::fabs(
s[
i] - s_node[
i]) > tol)
3962 unsigned nnod =
nnode();
3963 for (
unsigned j = 0; j < nnod; j++)
3972 cog[
i] /= double(nnod);
3976 for (
unsigned j = 0; j < nnod; j++)
3978 double dist_squared = 0.0;
3984 if (dist_squared > max_radius) max_radius = dist_squared;
3986 max_radius = sqrt(max_radius);
3993 const unsigned&
i)
const
3996 const unsigned n_node =
nnode();
4000 Shape psi(n_node, n_position_type);
4007 for (
unsigned l = 0; l < n_node; l++)
4010 for (
unsigned k = 0; k < n_position_type; k++)
4025 const unsigned&
i)
const
4028 const unsigned n_node =
nnode();
4033 Shape psi(n_node, n_position_type);
4040 for (
unsigned l = 0; l < n_node; l++)
4043 for (
unsigned k = 0; k < n_position_type; k++)
4059 const unsigned n_node =
nnode();
4066 Shape psi(n_node, n_position_type);
4071 for (
unsigned i = 0;
i < nodal_dim;
i++)
4076 for (
unsigned l = 0; l < n_node; l++)
4079 for (
unsigned k = 0; k < n_position_type; k++)
4096 const unsigned n_node =
nnode();
4103 Shape psi(n_node, n_position_type);
4108 for (
unsigned i = 0;
i < nodal_dim;
i++)
4113 for (
unsigned l = 0; l < n_node; l++)
4116 for (
unsigned k = 0; k < n_position_type; k++)
4136 const unsigned n_node =
nnode();
4140 const unsigned n_dim_element =
dim();
4143 Shape psi(n_node, n_position_type);
4144 DShape dpsids(n_node, n_position_type, n_dim_element);
4154 for (
unsigned i = 0;
i < n_dim_element;
i++)
4156 for (
unsigned j = 0; j < n_dim_element; j++)
4158 for (
unsigned k = 0; k < n_dim_node; k++)
4160 G(
i, j) += interpolated_G(
i, k) * interpolated_G(j, k);
4167 switch (n_dim_element)
4170 throw OomphLibError(
"Cannot calculate J_eulerian() for point element\n",
4171 OOMPH_CURRENT_FUNCTION,
4172 OOMPH_EXCEPTION_LOCATION);
4178 det = G(0, 0) * G(1, 1) - G(0, 1) * G(1, 0);
4181 det = G(0, 0) * G(1, 1) * G(2, 2) + G(0, 1) * G(1, 2) * G(2, 0) +
4182 G(0, 2) * G(1, 0) * G(2, 1) - G(0, 0) * G(1, 2) * G(2, 1) -
4183 G(0, 1) * G(1, 0) * G(2, 2) - G(0, 2) * G(1, 1) * G(2, 0);
4186 oomph_info <<
"More than 3 dimensions in J_eulerian()" << std::endl;
4201 const unsigned n_node =
nnode();
4206 const unsigned n_dim_element =
dim();
4209 Shape psi(n_node, n_position_type);
4210 DShape dpsids(n_node, n_position_type, n_dim_element);
4223 for (
unsigned i = 0;
i < n_dim_element;
i++)
4225 for (
unsigned j = 0; j < n_dim_element; j++)
4227 for (
unsigned k = 0; k < n_dim_node; k++)
4229 G(
i, j) += interpolated_G(
i, k) * interpolated_G(j, k);
4236 switch (n_dim_element)
4239 throw OomphLibError(
"Cannot calculate J_eulerian() for point element\n",
4240 OOMPH_CURRENT_FUNCTION,
4241 OOMPH_EXCEPTION_LOCATION);
4247 det = G(0, 0) * G(1, 1) - G(0, 1) * G(1, 0);
4250 det = G(0, 0) * G(1, 1) * G(2, 2) + G(0, 1) * G(1, 2) * G(2, 0) +
4251 G(0, 2) * G(1, 0) * G(2, 1) - G(0, 0) * G(1, 2) * G(2, 1) -
4252 G(0, 1) * G(1, 0) * G(2, 2) - G(0, 2) * G(1, 1) * G(2, 0);
4255 oomph_info <<
"More than 3 dimensions in J_eulerian()" << std::endl;
4277 const unsigned n_node =
nnode();
4286 const unsigned n_dim_element =
dim();
4289 Shape psi(n_node, n_position_type);
4290 DShape dpsi(n_node, n_dim_element);
4293 for (
unsigned ipt = 0; ipt < nintpt; ipt++)
4327 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
4351 const unsigned ncomponents = integral.size();
4353 for (
unsigned i = 0;
i < ncomponents;
i++)
4369 const unsigned n_dim =
dim();
4373 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
4376 for (
unsigned i = 0;
i < n_dim;
i++)
4382 for (
unsigned i = 0;
i < n_dim_eulerian;
i++)
4394 integrand_fct_pt(time, x, integrand);
4397 for (
unsigned i = 0;
i < ncomponents;
i++)
4399 integral[
i] += integrand[
i] * w * J;
4412 const unsigned ncomponents = integral.size();
4414 for (
unsigned i = 0;
i < ncomponents;
i++)
4430 const unsigned n_dim =
dim();
4434 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
4437 for (
unsigned i = 0;
i < n_dim;
i++)
4443 for (
unsigned i = 0;
i < n_dim_eulerian;
i++)
4455 integrand_fct_pt(x, integrand);
4458 for (
unsigned i = 0;
i < ncomponents;
i++)
4460 integral[
i] += integrand[
i] * w * J;
4485 OomphLibWarning(
"Pointer to spatial integration scheme has not been set.",
4486 "FiniteElement::self_test()",
4487 OOMPH_EXCEPTION_LOCATION);
4492 const unsigned dim_el =
dim();
4506 const unsigned n_node =
nnode();
4509 Shape psi(n_node, n_position_type);
4510 DShape dpsidx(n_node, dim_el);
4523 if (face_el_pt == 0)
4533 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
4536 for (
unsigned i = 0;
i < dim_el;
i++)
4543 for (
unsigned test = 0; test < ntest; test++)
4567 OOMPH_CURRENT_FUNCTION,
4568 OOMPH_EXCEPTION_LOCATION);
4573 if (std::fabs(jacobian) < 1.0e-16)
4575 std::ostringstream warning_stream;
4576 warning_stream <<
"Determinant of Jacobian matrix is zero at ipt "
4577 << ipt << std::endl;
4579 "FiniteElement::self_test()",
4580 OOMPH_EXCEPTION_LOCATION);
4589 std::ostringstream warning_stream;
4590 warning_stream <<
"Jacobian negative at integration point ipt="
4591 << ipt << std::endl;
4593 <<
"If you think that this is what you want you may: "
4595 <<
"set the (static) flag "
4596 <<
"FiniteElement::Accept_negative_jacobian to be true"
4600 "FiniteElement::self_test()",
4601 OOMPH_EXCEPTION_LOCATION);
4628 const unsigned& t_deriv)
4631 const unsigned n_node =
nnode();
4634 Shape psi(n_node, n_position_type);
4642 for (
unsigned l = 0; l < n_node; l++)
4645 for (
unsigned k = 0; k < n_position_type; k++)
4660 const unsigned& t_deriv,
4664 const unsigned n_node =
nnode();
4669 Shape psi(n_node, n_position_type);
4673 for (
unsigned i = 0;
i < nodal_dim;
i++)
4680 for (
unsigned l = 0; l < n_node; l++)
4683 for (
unsigned k = 0; k < n_position_type; k++)
4717 const unsigned n_node = this->
nnode();
4721 Shape psi(n_node, n_position_type);
4723 this->
shape(s, psi);
4726 const unsigned ncoord = this->
dim();
4728 for (
unsigned i = 0;
i < ncoord;
i++)
4735 for (
unsigned l = 0; l < n_node; l++)
4737 for (
unsigned k = 0; k < n_position_type; k++)
4740 const double psi_ = psi(l, k);
4741 for (
unsigned i = 0;
i < ncoord;
i++)
4767 const bool& use_coordinate_as_initial_guess)
4771 unsigned ncoord = this->
dim();
4779 double max_radius = 0.0;
4783 double radius = 0.0;
4784 for (
unsigned i = 0;
i < ncoord;
i++)
4786 radius += (cog[
i] - zeta[
i]) * (cog[
i] - zeta[
i]);
4788 radius = sqrt(radius);
4789 if (radius > Locate_zeta_helpers::
4812 if (!use_coordinate_as_initial_guess)
4818 unsigned num_plot_points =
4820 for (
unsigned iplot = 0; iplot < num_plot_points; iplot++)
4824 s_list.push_back(s_c);
4832 bool keep_going =
true;
4838 if (use_coordinate_as_initial_guess)
4844 for (
unsigned i = 0;
i < ncoord;
i++)
4846 dx[
i] = zeta[
i] - inter_x[
i];
4852 double my_min_resid = DBL_MAX;
4857 unsigned n_list_coord = s_list.size();
4859 for (
unsigned i_coord = 0; i_coord < n_list_coord; i_coord++)
4861 for (
unsigned i = 0;
i < ncoord;
i++)
4863 s_local[
i] = s_list[i_coord][
i];
4869 for (
unsigned i = 0;
i < ncoord;
i++)
4871 work_dx[
i] = zeta[
i] - work_x[
i];
4874 double maxres = std::fabs(
4875 *std::max_element(work_dx.begin(), work_dx.end(),
AbsCmp<double>()));
4878 if (maxres < my_min_resid)
4880 my_min_resid = maxres;
4905 std::fabs(*std::max_element(dx.begin(), dx.end(),
AbsCmp<double>()));
4929 for (
unsigned i = 0;
i < ncoord;
i++)
4932 work_s[
i] += fd_step;
4938 for (
unsigned j = 0; j < ncoord; j++)
4940 jacobian(j,
i) = -(work_r[j] - r[j]) / fd_step;
4950 unsigned n_node = this->
nnode();
4952 Shape psi(n_node, n_position_type);
4953 DShape dpsids(n_node, n_position_type, ncoord);
4982 for (
unsigned l = 0; l < n_node; l++)
4986 for (
unsigned k = 0; k < n_position_type; k++)
4989 for (
unsigned i = 0;
i < ncoord;
i++)
4991 for (
unsigned j = 0; j < ncoord; j++)
4993 interpolated_dxds(
i, j) +=
5002 for (
unsigned i = 0;
i < ncoord;
i++)
5004 for (
unsigned j = 0; j < ncoord; j++)
5006 jacobian(
i, j) = -interpolated_dxds(
i, j);
5022 <<
"FiniteElement::locate_zeta()" << std::endl;
5023 oomph_info <<
"Should not affect the result!" << std::endl;
5028 for (
unsigned i = 0;
i < ncoord;
i++)
5035 for (
unsigned i = 0;
i < ncoord;
i++)
5037 dx[
i] = zeta[
i] - inter_x[
i];
5042 std::fabs(*std::max_element(dx.begin(), dx.end(),
AbsCmp<double>()));
5050 }
while (keep_going);
5063 for (
unsigned i = 0;
i < ncoord;
i++)
5065 dx[
i] = zeta[
i] - inter_x[
i];
5070 std::fabs(*std::max_element(dx.begin(), dx.end(),
AbsCmp<double>()));
5094 geom_object_pt =
this;
5104 const unsigned n_node =
nnode();
5105 for (
unsigned n = 0; n < n_node; n++)
5127 std::set<std::pair<Data*, unsigned>>& paired_field_data)
5131 for (
unsigned n = 0; n < n_internal; n++)
5136 const unsigned n_value = dat_pt->
nvalue();
5139 for (
unsigned i = 0;
i < n_value;
i++)
5141 paired_field_data.insert(std::make_pair(dat_pt,
i));
5146 const unsigned n_node = this->
nnode();
5147 for (
unsigned n = 0; n < n_node; n++)
5152 const unsigned n_value = nod_pt->
nvalue();
5155 for (
unsigned i = 0;
i < n_value;
i++)
5157 paired_field_data.insert(std::make_pair(nod_pt,
i));
5171 #ifdef OOMPH_HAS_MPI
5201 for (
unsigned i = 0;
i < nnode_face;
i++)
5233 unsigned n_dim =
dim();
5241 for (
unsigned iplot = 0; iplot < num_plot_points; iplot++)
5247 for (
unsigned i = 0;
i < n_dim + 1;
i++)
5255 for (
unsigned i = 0;
i < n_dim;
i++)
5257 outfile << zeta[
i] <<
" ";
5259 outfile << std::endl;
5275 unsigned n_dim_el = this->
dim();
5279 if (n_dim_el == 0)
return 1.0;
5282 unsigned n_node =
nnode();
5291 Shape psi(n_node, n_position_type);
5292 DShape dpsids(n_node, n_position_type, n_dim_el);
5301 for (
unsigned l = 0; l < n_node; l++)
5304 for (
unsigned k = 0; k < n_position_type; k++)
5307 for (
unsigned i = 0;
i < n_dim;
i++)
5310 for (
unsigned j = 0; j < n_dim_el; j++)
5312 interpolated_A(j,
i) +=
5320 for (
unsigned i = 0;
i < n_dim_el;
i++)
5322 for (
unsigned j = 0; j < n_dim_el; j++)
5325 for (
unsigned k = 0; k < n_dim; k++)
5327 A(
i, j) += interpolated_A(
i, k) * interpolated_A(j, k);
5340 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
5344 "FaceElement::J_eulerian()",
5345 OOMPH_EXCEPTION_LOCATION);
5361 const unsigned n_node =
nnode();
5368 const unsigned n_dim_el =
dim();
5371 Shape psi(n_node, n_position_type);
5372 DShape dpsids(n_node, n_position_type, n_dim_el);
5384 for (
unsigned l = 0; l < n_node; l++)
5387 for (
unsigned k = 0; k < n_position_type; k++)
5390 for (
unsigned i = 0;
i < n_dim;
i++)
5393 for (
unsigned j = 0; j < n_dim_el; j++)
5395 interpolated_A(j,
i) +=
5404 for (
unsigned i = 0;
i < n_dim_el;
i++)
5406 for (
unsigned j = 0; j < n_dim_el; j++)
5409 for (
unsigned k = 0; k < n_dim; k++)
5411 A(
i, j) += interpolated_A(
i, k) * interpolated_A(j, k);
5424 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
5428 "FaceElement::J_eulerian_at_knot()",
5429 OOMPH_EXCEPTION_LOCATION);
5446 const unsigned n_node =
nnode();
5453 const unsigned n_dim_el =
dim();
5456 Shape psi(n_node, n_position_type);
5457 DShape dpsids(n_node, n_position_type, n_dim_el);
5461 for (
unsigned ipt = 0; ipt < nintpt; ipt++)
5473 for (
unsigned l = 0; l < n_node; l++)
5476 for (
unsigned k = 0; k < n_position_type; k++)
5479 for (
unsigned i = 0;
i < n_dim;
i++)
5483 for (
unsigned j = 0; j < n_dim_el; j++)
5485 interpolated_A(j,
i) +=
5495 for (
unsigned i = 0;
i < n_dim_el;
i++)
5497 for (
unsigned j = 0; j < n_dim_el; j++)
5500 for (
unsigned k = 0; k < n_dim; k++)
5502 A(
i, j) += interpolated_A(
i, k) * interpolated_A(j, k);
5515 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
5519 "FaceElement::J_eulerian_at_knot()",
5520 OOMPH_EXCEPTION_LOCATION);
5548 const unsigned element_dim =
dim();
5556 if (
s.size() != element_dim)
5558 std::ostringstream error_stream;
5560 <<
"Local coordinate s passed to outer_unit_normal() has dimension "
5561 <<
s.size() << std::endl
5562 <<
"but element dimension is " << element_dim << std::endl;
5565 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5573 std::ostringstream error_stream;
5574 error_stream <<
"Tangent direction vector has dimension "
5576 <<
"but spatial dimension is " << spatial_dim << std::endl;
5579 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5583 if (unit_normal.size() != spatial_dim)
5585 std::ostringstream error_stream;
5586 error_stream <<
"Unit normal passed to outer_unit_normal() has dimension "
5587 << unit_normal.size() << std::endl
5588 <<
"but spatial dimension is " << spatial_dim << std::endl;
5591 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5596 unsigned ntangent_vec = tang_vec.size();
5602 switch (element_dim)
5608 if (ntangent_vec != 1)
5610 std::ostringstream error_stream;
5612 <<
"This is a 0D FaceElement, we need one tangent vector.\n"
5613 <<
"You have given me " << tang_vec.size() <<
" vector(s).\n";
5615 OOMPH_CURRENT_FUNCTION,
5616 OOMPH_EXCEPTION_LOCATION);
5625 if (ntangent_vec != 1)
5627 std::ostringstream error_stream;
5629 <<
"This is a 1D FaceElement, we need one tangent vector.\n"
5630 <<
"You have given me " << tang_vec.size() <<
" vector(s).\n";
5632 OOMPH_CURRENT_FUNCTION,
5633 OOMPH_EXCEPTION_LOCATION);
5642 if (ntangent_vec != 2)
5644 std::ostringstream error_stream;
5646 <<
"This is a 2D FaceElement, we need two tangent vectors.\n"
5647 <<
"You have given me " << tang_vec.size() <<
" vector(s).\n";
5649 OOMPH_CURRENT_FUNCTION,
5650 OOMPH_EXCEPTION_LOCATION);
5658 "Cannot have a FaceElement with dimension higher than 2",
5659 OOMPH_CURRENT_FUNCTION,
5660 OOMPH_EXCEPTION_LOCATION);
5665 for (
unsigned vec_i = 0; vec_i < ntangent_vec; vec_i++)
5667 if (tang_vec[vec_i].
size() != spatial_dim)
5669 std::ostringstream error_stream;
5670 error_stream <<
"This problem has " << spatial_dim
5671 <<
" spatial dimensions.\n"
5672 <<
"But the " << vec_i <<
" tangent vector has size "
5673 << tang_vec[vec_i].size() << std::endl;
5675 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5683 switch (element_dim)
5688 std::ostringstream error_stream;
5690 <<
"It is unclear how to calculate a tangent and normal vectors for "
5691 <<
"a point element.\n";
5693 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5706 const unsigned n_position_type_bulk =
5716 Shape psi(n_node_bulk, n_position_type_bulk);
5717 DShape dpsids(n_node_bulk, n_position_type_bulk, 2);
5726 for (
unsigned l = 0; l < n_node_bulk; l++)
5729 for (
unsigned k = 0; k < n_position_type_bulk; k++)
5732 for (
unsigned j = 0; j < 2; j++)
5735 for (
unsigned i = 0;
i < spatial_dim;
i++)
5738 interpolated_dxds(j,
i) +=
5758 unsigned interior_direction = 0;
5761 for (
unsigned i = 0;
i < spatial_dim;
i++)
5765 tangent[
i] = interpolated_dxds(0,
i) * dsbulk_dsface(0, 0) +
5766 interpolated_dxds(1,
i) * dsbulk_dsface(1, 0);
5769 interior_tangent[
i] = interpolated_dxds(interior_direction,
i);
5783 for (
unsigned i = 0;
i < spatial_dim;
i++)
5785 unit_normal[
i] = full_normal[
i];
5790 for (
unsigned i = 0;
i < spatial_dim;
i++)
5796 tang_vec[0][0] = -unit_normal[1];
5797 tang_vec[0][1] = unit_normal[0];
5811 if (spatial_dim != 3)
5813 std::ostringstream error_stream;
5814 error_stream <<
"There are only " << spatial_dim
5815 <<
"coordinates at the nodes of this 2D FaceElement,\n"
5816 <<
"which must have come from a 3D Bulk element\n";
5818 OOMPH_CURRENT_FUNCTION,
5819 OOMPH_EXCEPTION_LOCATION);
5824 const unsigned n_node = this->
nnode();
5830 Shape psi(n_node, n_position_type);
5831 DShape dpsids(n_node, n_position_type, 2);
5840 for (
unsigned l = 0; l < n_node; l++)
5843 for (
unsigned k = 0; k < n_position_type; k++)
5846 for (
unsigned j = 0; j < 2; j++)
5849 for (
unsigned i = 0;
i < 3;
i++)
5854 interpolated_dxds[j][
i] +=
5864 interpolated_dxds[0], interpolated_dxds[1], unit_normal);
5869 for (
unsigned i = 0;
i < spatial_dim;
i++)
5894 std::ostringstream err_stream;
5895 err_stream <<
"The angle between the unit normal and the \n"
5896 <<
"direction vector is less than 10 degrees.";
5898 OOMPH_CURRENT_FUNCTION,
5899 OOMPH_EXCEPTION_LOCATION);
5913 tang_vec[0][0] = t1_scaling * interpolated_dxds[0][0] +
5914 t2_scaling * interpolated_dxds[1][0];
5915 tang_vec[0][1] = t1_scaling * interpolated_dxds[0][1] +
5916 t2_scaling * interpolated_dxds[1][1];
5917 tang_vec[0][2] = t1_scaling * interpolated_dxds[0][2] +
5918 t2_scaling * interpolated_dxds[1][2];
5925 for (
unsigned vec_i = 0; vec_i < 2; vec_i++)
5930 for (
unsigned dim_i = 0; dim_i < spatial_dim; dim_i++)
5932 tang_vec[vec_i][dim_i] /= tang_length;
5949 if (a != 0.0 || b != 0.0)
5951 double a_sq = a * a;
5952 double b_sq = b * b;
5953 double c_sq = c * c;
5955 tang_vec[0][0] = -b / sqrt(a_sq + b_sq);
5956 tang_vec[0][1] = a / sqrt(a_sq + b_sq);
5959 double z = (a_sq + b_sq) / sqrt(a_sq * c_sq + b_sq * c_sq +
5960 (a_sq + b_sq) * (a_sq + b_sq));
5962 tang_vec[1][0] = -(a * c * z) / (a_sq + b_sq);
5963 tang_vec[1][1] = -(b * c * z) / (a_sq + b_sq);
5972 std::ostringstream warn_stream;
5974 <<
"I have detected that your normal vector is [0,0,1].\n"
5975 <<
"Since this function compares against 0.0, you may\n"
5976 <<
"get discontinuous tangent vectors.";
5978 OOMPH_CURRENT_FUNCTION,
5979 OOMPH_EXCEPTION_LOCATION);
5982 tang_vec[0][0] = 1.0;
5983 tang_vec[0][1] = 0.0;
5984 tang_vec[0][2] = 0.0;
5986 tang_vec[1][0] = 0.0;
5987 tang_vec[1][1] = 1.0;
5988 tang_vec[1][2] = 0.0;
5993 OOMPH_CURRENT_FUNCTION,
5994 OOMPH_EXCEPTION_LOCATION);
6003 "Cannot have a FaceElement with dimension higher than 2",
6004 OOMPH_CURRENT_FUNCTION,
6005 OOMPH_EXCEPTION_LOCATION);
6017 const unsigned& ipt,
6022 const unsigned element_dim =
dim();
6025 for (
unsigned i = 0;
i < element_dim;
i++)
6040 const unsigned element_dim =
dim();
6048 if (
s.size() != element_dim)
6050 std::ostringstream error_stream;
6052 <<
"Local coordinate s passed to outer_unit_normal() has dimension "
6053 <<
s.size() << std::endl
6054 <<
"but element dimension is " << element_dim << std::endl;
6057 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
6061 if (unit_normal.size() != spatial_dim)
6063 std::ostringstream error_stream;
6064 error_stream <<
"Unit normal passed to outer_unit_normal() has dimension "
6065 << unit_normal.size() << std::endl
6066 <<
"but spatial dimension is " << spatial_dim << std::endl;
6069 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
6126 switch (element_dim)
6138 const unsigned n_position_type_bulk =
6149 Shape psi(n_node_bulk, n_position_type_bulk);
6150 DShape dpsids(n_node_bulk, n_position_type_bulk, 1);
6159 for (
unsigned l = 0; l < n_node_bulk; l++)
6162 for (
unsigned k = 0; k < n_position_type_bulk; k++)
6165 for (
unsigned i = 0;
i < spatial_dim;
i++)
6168 interpolated_dxds(0,
i) +=
6176 for (
unsigned i = 0;
i < spatial_dim;
i++)
6178 unit_normal[
i] = interpolated_dxds(0,
i);
6192 const unsigned n_position_type_bulk =
6202 Shape psi(n_node_bulk, n_position_type_bulk);
6203 DShape dpsids(n_node_bulk, n_position_type_bulk, 2);
6212 for (
unsigned l = 0; l < n_node_bulk; l++)
6215 for (
unsigned k = 0; k < n_position_type_bulk; k++)
6218 for (
unsigned j = 0; j < 2; j++)
6221 for (
unsigned i = 0;
i < spatial_dim;
i++)
6224 interpolated_dxds(j,
i) +=
6244 unsigned interior_direction = 0;
6247 for (
unsigned i = 0;
i < spatial_dim;
i++)
6251 tangent[
i] = interpolated_dxds(0,
i) * dsbulk_dsface(0, 0) +
6252 interpolated_dxds(1,
i) * dsbulk_dsface(1, 0);
6255 interior_tangent[
i] = interpolated_dxds(interior_direction,
i);
6269 for (
unsigned i = 0;
i < spatial_dim;
i++)
6271 unit_normal[
i] = full_normal[
i];
6286 if (spatial_dim != 3)
6288 std::ostringstream error_stream;
6289 error_stream <<
"There are only " << spatial_dim
6290 <<
"coordinates at the nodes of this 2D FaceElement,\n"
6291 <<
"which must have come from a 3D Bulk element\n";
6293 OOMPH_CURRENT_FUNCTION,
6294 OOMPH_EXCEPTION_LOCATION);
6299 const unsigned n_node = this->
nnode();
6305 Shape psi(n_node, n_position_type);
6306 DShape dpsids(n_node, n_position_type, 2);
6315 for (
unsigned l = 0; l < n_node; l++)
6318 for (
unsigned k = 0; k < n_position_type; k++)
6321 for (
unsigned j = 0; j < 2; j++)
6324 for (
unsigned i = 0;
i < 3;
i++)
6329 interpolated_dxds[j][
i] +=
6339 interpolated_dxds[0], interpolated_dxds[1], unit_normal);
6346 "Cannot have a FaceElement with dimension higher than 2",
6347 OOMPH_CURRENT_FUNCTION,
6348 OOMPH_EXCEPTION_LOCATION);
6354 for (
unsigned i = 0;
i < spatial_dim;
i++)
6367 const unsigned element_dim =
dim();
6370 for (
unsigned i = 0;
i < element_dim;
i++)
6396 (*Face_to_bulk_coordinate_fct_pt)(
s, s_bulk);
6400 throw OomphLibError(
"Face_to_bulk_coordinate mapping not set",
6401 OOMPH_CURRENT_FUNCTION,
6402 OOMPH_EXCEPTION_LOCATION);
6421 (*Face_to_bulk_coordinate_fct_pt)(
s, s_bulk);
6426 throw OomphLibError(
"Face_to_bulk_coordinate mapping not set",
6427 OOMPH_CURRENT_FUNCTION,
6428 OOMPH_EXCEPTION_LOCATION);
6441 unsigned& interior_direction)
const
6447 (*Bulk_coordinate_derivatives_fct_pt)(
6448 s, dsbulk_dsface, interior_direction);
6454 "No function for derivatives of bulk coords wrt face coords set",
6455 OOMPH_CURRENT_FUNCTION,
6456 OOMPH_EXCEPTION_LOCATION);
6477 unsigned n_dim =
dim();
6489 unsigned n_node = this->
nnode();
6498 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
6501 for (
unsigned i = 0;
i < n_dim;
i++)
6520 for (
unsigned idim = 0; idim < n_dim; idim++)
6522 disp[idim] = x[idim] - xi[idim];
6526 for (
unsigned ii = 0; ii < n_dim; ii++)
6528 el_norm += (disp[ii] * disp[ii]) *
W;
6545 std::ostream& out,
const std::string& current_string)
const
6562 const unsigned el_dim =
dim();
6566 const unsigned n_shape =
nnode();
6573 std::ostringstream error_message;
6574 error_message <<
"Dimension mismatch" << std::endl;
6575 error_message <<
"The elemental dimension: " << el_dim
6576 <<
" must equal the nodal Lagrangian dimension: "
6578 <<
" for the jacobian of the mapping to be well-defined"
6581 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
6586 for (
unsigned i = 0;
i < el_dim;
i++)
6589 for (
unsigned j = 0; j < el_dim; j++)
6592 jacobian(
i, j) = 0.0;
6594 for (
unsigned l = 0; l < n_shape; l++)
6596 for (
unsigned k = 0; k < n_shape_type; k++)
6620 const unsigned el_dim =
dim();
6624 const unsigned n_shape =
nnode();
6627 const unsigned n_row =
N2deriv[el_dim];
6632 for (
unsigned i = 0;
i < n_row;
i++)
6635 for (
unsigned j = 0; j < el_dim; j++)
6638 jacobian2(
i, j) = 0.0;
6640 for (
unsigned l = 0; l < n_shape; l++)
6643 for (
unsigned k = 0; k < n_shape_type; k++)
6679 const unsigned el_dim =
dim();
6683 const unsigned n_shape =
nnode();
6690 for (
unsigned i = 0;
i < el_dim;
i++)
6693 for (
unsigned j = 0; j < el_dim; j++)
6695 jacobian(
i, j) = 0.0;
6696 inverse_jacobian(
i, j) = 0.0;
6700 for (
unsigned l = 0; l < n_shape; l++)
6703 for (
unsigned k = 0; k < n_shape_type; k++)
6714 for (
unsigned i = 0;
i < el_dim;
i++)
6716 det *= jacobian(
i,
i);
6725 for (
unsigned i = 0;
i < el_dim;
i++)
6727 inverse_jacobian(
i,
i) = 1.0 / jacobian(
i,
i);
6745 const unsigned el_dim =
dim();
6772 const unsigned el_dim =
dim();
6817 const unsigned el_dim =
dim();
6819 const unsigned n_deriv =
N2deriv[el_dim];
6837 jacobian, inverse_jacobian, jacobian2, dpsi, d2psi);
6869 const unsigned el_dim =
dim();
6871 const unsigned n_deriv =
N2deriv[el_dim];
6889 jacobian, inverse_jacobian, jacobian2, dpsi, d2psi);
6905 std::ostream& out,
const std::string& current_string)
const
6908 const unsigned n_node =
nnode();
6910 for (
unsigned n = 0; n < n_node; n++)
6914 std::stringstream conversion;
6915 conversion <<
" of Solid Node " << n << current_string;
6929 const bool& store_local_dof_pt)
6932 const unsigned n_node =
nnode();
6948 std::deque<unsigned long> global_eqn_number_queue;
6955 for (
unsigned n = 0; n < n_node; n++)
6961 for (
unsigned j = 0; j < n_position_type; j++)
6964 for (
unsigned k = 0; k < nodal_dim; k++)
6973 global_eqn_number_queue.push_back(
eqn_number);
6975 if (store_local_dof_pt)
6978 &(cast_node_pt->
x_gen(j, k)));
7000 if (store_local_dof_pt)
7022 const unsigned n_node =
nnode();
7039 const unsigned n_dof = this->
ndof();
7046 int local_unknown = 0;
7052 for (
unsigned n = 0; n < n_node; n++)
7055 for (
unsigned k = 0; k < n_position_type; k++)
7058 for (
unsigned i = 0;
i < nodal_dim;
i++)
7062 if (local_unknown >= 0)
7068 const double old_var = *value_pt;
7071 *value_pt += fd_step;
7085 for (
unsigned m = 0; m < n_dof; m++)
7088 jacobian(m, local_unknown) =
7089 (newres[m] - residuals[m]) / fd_step;
7112 *value_pt = old_var;
7135 const unsigned&
i)
const
7138 const unsigned n_node =
nnode();
7144 Shape psi(n_node, n_lagrangian_type);
7153 for (
unsigned l = 0; l < n_node; l++)
7156 for (
unsigned k = 0; k < n_lagrangian_type; k++)
7173 const unsigned n_node =
nnode();
7179 Shape psi(n_node, n_lagrangian_type);
7188 for (
unsigned i = 0;
i < n_lagrangian;
i++)
7194 for (
unsigned l = 0; l < n_node; l++)
7197 for (
unsigned k = 0; k < n_lagrangian_type; k++)
7214 const unsigned n_node =
nnode();
7220 unsigned el_dim =
dim();
7223 Shape psi(n_node, n_lagrangian_type);
7224 DShape dpsi(n_node, n_lagrangian_type, el_dim);
7233 for (
unsigned i = 0;
i < n_lagrangian;
i++)
7235 for (
unsigned j = 0; j < el_dim; j++)
7241 for (
unsigned l = 0; l < n_node; l++)
7244 for (
unsigned k = 0; k < n_lagrangian_type; k++)
7265 std::ostringstream error_stream;
7266 error_stream <<
"Can only call fill_in_jacobian_for_newmark_accel()\n"
7267 <<
"With Solve_for_consistent_newmark_accel_flag:"
7269 error_stream <<
"Solid_ic_pt: " <<
Solid_ic_pt << std::endl;
7272 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
7277 const unsigned n_node =
nnode();
7282 const unsigned n_lagrangian =
dim();
7285 int local_eqn = 0, local_unknown = 0;
7289 Shape psi(n_node, n_position_type);
7293 DShape dpsidxi(n_node, n_position_type, n_lagrangian);
7307 if (timestepper_pt->
type() !=
"Newmark")
7309 std::ostringstream error_stream;
7311 <<
"Assignment of Newmark accelerations obviously only works\n"
7312 <<
"for Newmark timesteppers. You've called me with "
7313 << timestepper_pt->
type() << std::endl;
7316 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
7321 const unsigned ntstorage = timestepper_pt->
ntstorage();
7322 const double accel_weight = timestepper_pt->
weight(2, ntstorage - 1);
7325 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
7328 for (
unsigned i = 0;
i < n_lagrangian;
i++)
7351 for (
unsigned l0 = 0; l0 < n_node; l0++)
7354 for (
unsigned k0 = 0; k0 < n_position_type; k0++)
7357 for (
unsigned i0 = 0; i0 < nodal_dim; i0++)
7364 for (
unsigned l1 = 0; l1 < n_node; l1++)
7368 for (
unsigned k1 = 0; k1 < n_position_type; k1++)
7375 if (local_unknown >= 0)
7380 jacobian(local_eqn, local_unknown) +=
7381 factor * accel_weight * psi(l0, k0) * psi(l1, k1) *
W;
7409 const unsigned& flag)
7412 const unsigned n_node =
nnode();
7422 const unsigned n_lagrangian =
dim();
7426 int local_unknown = 0;
7429 Shape psi(n_node, n_position_type);
7442 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
7445 for (
unsigned i = 0;
i < n_lagrangian;
i++)
7460 for (
unsigned i = 0;
i < n_lagrangian;
i++)
7463 for (
unsigned l = 0; l < n_node; l++)
7466 for (
unsigned k = 0; k < n_lagrangian_type; k++)
7481 for (
unsigned l = 0; l < n_node; l++)
7484 for (
unsigned k = 0; k < n_position_type; k++)
7487 for (
unsigned i = 0;
i < nodal_dim;
i++)
7498 residuals[local_eqn] +=
7506 for (
unsigned ll = 0; ll < n_node; ll++)
7509 for (
unsigned kk = 0; kk < n_position_type; kk++)
7517 if (local_unknown >= 0)
7523 jacobian(local_eqn, local_unknown) +=
7524 psi(ll, kk) * psi(l, k) * w;
7529 <<
"WARNING: You should really free all Data"
7531 <<
" before setup of initial guess" << std::endl
7532 <<
"ll, kk, ii " << ll <<
" " << kk <<
" " << ii
7541 oomph_info <<
"WARNING: You should really free all Data"
7543 <<
" before setup of initial guess"
7545 <<
"l, k, i " << l <<
" " << k <<
" " <<
i
Function-type-object to perform absolute comparison of objects. Apparently this inlines better.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
unsigned long nindex2() const
Return the range of index 2 of the derivatives of the shape functions.
unsigned long nindex1() const
Return the range of index 1 of the derivatives of the shape functions.
A class that represents a collection of data; each Data object may contain many different individual ...
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
virtual void add_eqn_numbers_to_vector(Vector< long > &vector_of_eqn_numbers)
Add all equation numbers to the vector in the internal storage order.
virtual void describe_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the dofs of the Node. The ostream specifies the output stream to which the descr...
virtual void add_value_pt_to_map(std::map< unsigned, double * > &map_of_value_pt)
Add pointers to all unpinned and unconstrained data to a map indexed by (global) equation number.
static long Is_pinned
Static "Magic number" used in place of the equation number to indicate that the value is pinned.
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 read_values_from_vector(const Vector< double > &vector_of_values, unsigned &index)
Read all data and time history values from the vector starting from index. On return the index will b...
virtual void add_values_to_vector(Vector< double > &vector_of_values)
Add all data and time history values to the vector in the internal storage order.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
virtual void assign_eqn_numbers(unsigned long &global_ndof, Vector< double * > &dof_pt)
Assign global equation numbers; increment global number of unknowns, global_ndof; and add any new dof...
static long Is_unclassified
Static "Magic number" used in place of the equation number to denote a value that hasn't been classif...
virtual void read_eqn_numbers_from_vector(const Vector< long > &vector_of_eqn_numbers, unsigned &index)
Read all equation numbers from the vector starting from index. On return the index will be set to the...
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
unsigned long nrow() const
Return the number of rows of the matrix.
unsigned long ncol() const
Return the number of columns of the matrix.
void solve(DoubleVector &rhs)
Complete LU solve (replaces matrix by its LU decomposition and overwrites RHS with solution)....
This is a base class for all elements that require external sources (e.g. FSI, multi-domain problems ...
unsigned nexternal_interaction_geometric_data() const
Return the number of geometric Data items that affect the external interactions in this element: i....
unsigned nexternal_interaction_field_data() const
Return the number of Data items that affect the external interactions in this element....
Vector< Data * > external_interaction_field_data_pt() const
Return vector of pointers to the field Data objects that affect the interactions on the element.
Vector< Data * > external_interaction_geometric_data_pt() const
Return vector of pointers to the geometric Data objects that affect the interactions on the element.
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
unsigned & bulk_position_type(const unsigned &i)
Return the position type in the "bulk" element that corresponds to position type i on the FaceElement...
CoordinateMappingFctPt & face_to_bulk_coordinate_fct_pt()
Return the pointer to the function that maps the face coordinate to the bulk coordinate.
FiniteElement * Bulk_element_pt
Pointer to the associated higher-dimensional "bulk" element.
void bulk_node_number_resize(const unsigned &i)
Resize the storage for the bulk node numbers.
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
unsigned & nbulk_value(const unsigned &n)
Return the number of values originally stored at local node n (before the FaceElement added additiona...
int Normal_sign
Sign of outer unit normal (relative to cross-products of tangent vectors in the corresponding "bulk" ...
BulkCoordinateDerivativesFctPt Bulk_coordinate_derivatives_fct_pt
Pointer to a function that returns the derivatives of the local "bulk" coordinates with respect to th...
Vector< double > local_coordinate_in_bulk(const Vector< double > &s) const
Return vector of local coordinates in bulk element, given the local coordinates in this FaceElement.
Vector< double > * Tangent_direction_pt
A general direction pointer for the tangent vectors. This is used in the function continuous_tangent_...
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s....
static bool Ignore_discontinuous_tangent_warning
Ignores the warning when the tangent vectors from continuous_tangent_and_outer_unit_normal may not be...
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
void check_J_eulerian_at_knots(bool &passed) const
Check that Jacobian of mapping between local and Eulerian coordinates at all integration points is po...
BulkCoordinateDerivativesFctPt & bulk_coordinate_derivatives_fct_pt()
Return the pointer to the function that returns the derivatives of the bulk coordinates wrt the face ...
void get_ds_bulk_ds_face(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction) const
Calculate the derivatives of the local coordinates in the bulk element with respect to the local coor...
void continuous_tangent_and_outer_unit_normal(const Vector< double > &s, Vector< Vector< double >> &tang_vec, Vector< double > &unit_normal) const
Compute the tangent vector(s) and the outer unit normal vector at the specified local coordinate....
void get_local_coordinate_in_bulk(const Vector< double > &s, Vector< double > &s_bulk) const
Calculate the vector of local coordinate in the bulk element given the local coordinates in this Face...
CoordinateMappingFctPt Face_to_bulk_coordinate_fct_pt
Pointer to a function that translates the face coordinate to the coordinate in the bulk element.
int & normal_sign()
Sign of outer unit normal (relative to cross-products of tangent vectors in the corresponding "bulk" ...
void output_zeta(std::ostream &outfile, const unsigned &nplot)
Output boundary coordinate zeta.
unsigned & bulk_node_number(const unsigned &n)
Return the bulk node number that corresponds to the n-th local node number.
void nbulk_value_resize(const unsigned &i)
Resize the storage for the number of values originally stored at the local nodes to i entries.
A general Finite Element class.
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s.
virtual int face_outer_unit_normal_sign(const int &face_index) const
Get the sign of the outer unit normal on the face given by face_index.
double d2shape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
Compute the geometric shape functions and also first and second derivatives w.r.t....
virtual void update_before_nodal_fd()
Function that is called before the finite differencing of any nodal data. This may be overloaded to u...
double dJ_eulerian_at_knot(const unsigned &ipt, Shape &psi, DenseMatrix< double > &djacobian_dX) const
Compute the geometric shape functions (psi) at integration point ipt. Return the determinant of the j...
void set_nodal_dimension(const unsigned &nodal_dim)
Set the dimension of the nodes in the element. This will typically only be required when constructing...
virtual void dshape_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsids) const
Return the geometric shape function and its derivative w.r.t. the local coordinates at the ipt-th int...
virtual void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size (broken virtual)
Integral * Integral_pt
Pointer to the spatial integration scheme.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
virtual void transform_derivatives(const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
Convert derivative w.r.t.local coordinates to derivatives w.r.t the coordinates used to assemble the ...
virtual unsigned get_bulk_node_number(const int &face_index, const unsigned &i) const
Get the number of the ith node on face face_index (in the bulk node vector).
virtual Node * get_node_at_local_coordinate(const Vector< double > &s) const
If there is a node at this local coordinate, return the pointer to the node.
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
virtual double s_min() const
Min value of local coordinate.
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
virtual void assign_nodal_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers for Data stored at the nodes Virtual so that it can be overloaded b...
virtual void local_fraction_of_node(const unsigned &j, Vector< double > &s_fraction)
Get the local fraction of the node j in the element A dumb, but correct default implementation is pro...
virtual double invert_jacobian_mapping(const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
A template-free interface that takes the matrix passed as jacobian and return its inverse in inverse_...
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
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...
double dnodal_position_gen_dt(const unsigned &n, const unsigned &k, const unsigned &i) const
i-th component of time derivative (velocity) of the generalised position, dx(k,i)/dt at local node n....
virtual 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...
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
static double Tolerance_for_singular_jacobian
Tolerance below which the jacobian is considered singular.
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 locate_zeta(const Vector< double > &zeta, GeomObject *&geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
For a given value of zeta, the "global" intrinsic coordinate of a mesh of FiniteElements represented ...
void check_jacobian(const double &jacobian) const
Helper function used to check for singular or negative Jacobians in the transform from local to globa...
virtual void d2shape_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsids, DShape &d2psids) const
Return the geometric shape function and its first and second derivatives w.r.t. the local coordinates...
virtual void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Function to compute the geometric shape functions and also first and second derivatives w....
virtual unsigned required_nvalue(const unsigned &n) const
Number of values that must be stored at local node n by the element. The default is 0,...
static const double Node_location_tolerance
Default value for the tolerance to be used when locating nodes via local coordinates.
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
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...
MacroElement * Macro_elem_pt
Pointer to the element's macro element (NULL by default)
virtual 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...
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 double s_max() const
Max. value of local coordinate.
unsigned Nodal_dimension
The spatial dimension of the nodes in the element. We assume that nodes have the same spatial dimensi...
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
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...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
virtual double interpolated_dxdt(const Vector< double > &s, const unsigned &i, const unsigned &t)
Return t-th time-derivative of the i-th FE-interpolated Eulerian coordinate at local coordinate s.
virtual 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...
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
virtual double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the values of the "global" intrinsic coordinate, zeta, of a compound geometric object (a mesh...
Node ** Node_pt
Storage for pointers to the nodes in the element.
virtual double local_to_eulerian_mapping(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 function...
virtual BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt(const int &face_index) const
Get a pointer to the derivative of the mapping from face to bulk coordinates.
void get_centre_of_gravity_and_max_radius_in_terms_of_zeta(Vector< double > &cog, double &max_radius) const
Compute centre of gravity of all nodes and radius of node that is furthest from it....
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
virtual void reset_after_nodal_fd()
Function that is call after the finite differencing of the nodal data. This may be overloaded to rese...
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
unsigned Elemental_dimension
The spatial dimension of the element, i.e. the number of local coordinates used to parametrize it.
virtual bool local_coord_is_valid(const Vector< double > &s)
Broken assignment operator.
void check_J_eulerian_at_knots(bool &passed) const
Check that Jacobian of mapping between local and Eulerian coordinates at all integration points is po...
virtual void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
virtual void move_local_coord_back_into_element(Vector< double > &s) const
Adjust local coordinates so that they're located inside the element.
virtual 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...
virtual 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...
virtual void d_dshape_eulerian_dnodal_coordinates(const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
A template-free interface that calculates the derivative w.r.t. the nodal coordinates of the derivat...
virtual double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point.
MacroElement * macro_elem_pt()
Access function to pointer to macro element.
void integrate_fct(FiniteElement::SteadyExactSolutionFctPt integrand_fct_pt, Vector< double > &integral)
Evaluate integral of a Vector-valued function over the element.
virtual void describe_local_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the local dofs of the element[s]. The ostream specifies the output stream to whi...
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement....
virtual void transform_second_derivatives(const DenseMatrix< double > &jacobian, const DenseMatrix< double > &inverse_jacobian, const DenseMatrix< double > &jacobian2, DShape &dbasis, DShape &d2basis) const
Convert derivatives and second derivatives w.r.t. local coordiantes to derivatives and second derivat...
virtual void describe_nodal_local_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the local dofs of the element[s]. The ostream specifies the output stream to whi...
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
int ** Nodal_local_eqn
Storage for the local equation numbers associated with the values stored at the nodes.
virtual ~FiniteElement()
The destructor cleans up the static memory allocated for shape function storage. Internal and externa...
virtual void get_x_from_macro_element(const Vector< double > &s, Vector< double > &x) const
Global coordinates as function of local coordinates using macro element representation....
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation....
virtual void node_update()
Update the positions of all nodes in the element using each node update function. The default impleme...
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
unsigned Nnode
Number of nodes in the element.
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
static const unsigned N2deriv[]
Static array that holds the number of second derivatives as a function of the dimension of the elemen...
static bool Accept_negative_jacobian
Boolean that if set to true allows a negative jacobian in the transform between global and local coor...
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
virtual unsigned nnode_on_face() const
void transform_derivatives_diagonal(const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
Convert derivative w.r.t local coordinates to derivatives w.r.t the coordinates used to assemble the ...
double raw_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....
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
static bool Suppress_output_while_checking_for_inverted_elements
Static boolean to suppress output while checking for inverted elements.
int get_node_number(Node *const &node_pt) const
Return the number of the node *node_pt if this node is in the element, else return -1;.
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
virtual double d2shape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx, DShape &d2psidx) const
Return the geometric shape functions and also first and second derivatives w.r.t. global coordinates ...
virtual void dJ_eulerian_dnodal_coordinates(const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
A template-free interface that calculates the derivative of the jacobian of a mapping with respect to...
virtual CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt(const int &face_index) const
Get a pointer to the function mapping face coordinates to bulk coordinates.
static const unsigned Default_Initial_Nvalue
Default return value for required_nvalue(n) which gives the number of "data" values required by the e...
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK.
void fill_in_jacobian_from_external_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
Calculate the contributions to the jacobian from the external degrees of freedom using finite differe...
unsigned nexternal_data() const
Return the number of external data objects.
virtual void reset_in_external_fd(const unsigned &i)
Function called within the finite difference loop for external data after the values in the i-th exte...
virtual void fill_in_contribution_to_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
Add the elemental contribution to the mass matrix matrix. and the residuals vector....
virtual void update_in_internal_fd(const unsigned &i)
Function called within the finite difference loop for internal data after a change in any values in t...
bool is_halo() const
Is this element a halo?
virtual void reset_after_external_fd()
Function that is call after the finite differencing of the external data. This may be overloaded to r...
virtual void update_before_external_fd()
Function that is called before the finite differencing of any external data. This may be overloaded t...
void assign_internal_eqn_numbers(unsigned long &global_number, Vector< double * > &Dof_pt)
Assign the global equation numbers to the internal Data. The arguments are the current highest global...
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
static bool Suppress_warning_about_repeated_external_data
Static boolean to suppress warnings about repeated external data. Defaults to true.
bool internal_data_fd(const unsigned &i) const
Return the status of the boolean flag indicating whether the internal data is included in the finite ...
unsigned long * Eqn_number
Storage for the global equation numbers represented by the degrees of freedom in the element.
void read_internal_eqn_numbers_from_vector(const Vector< long > &vector_of_eqn_numbers, unsigned &index)
Read all equation numbers associated with internal data from the vector starting from index....
virtual void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Add the elemental contribution to the jacobian matrix, mass matrix and the residuals vector....
virtual void fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter(double *const ¶meter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
Add the elemental contribution to the derivative of the jacobian matrix, mass matrix and the residual...
virtual void fill_in_contribution_to_inner_product_vectors(Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
Fill in the contributions to the vectors that when taken as dot product with other history values giv...
void fill_in_jacobian_from_internal_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
Calculate the contributions to the jacobian from the internal degrees of freedom using finite differe...
bool external_data_fd(const unsigned &i) const
Return the status of the boolean flag indicating whether the external data is included in the finite ...
virtual void assign_all_generic_local_eqn_numbers(const bool &store_local_dof_pt)
Assign all the local equation numbering schemes that can be applied generically for the element....
unsigned add_internal_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) internal data object to the element and return the index required to obtain it ...
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.
int ** Data_local_eqn
Pointer to array storage for local equation numbers associated with internal and external data....
int local_eqn_number(const unsigned long &ieqn_global) const
Return the local equation number corresponding to the ieqn_global-th global equation number....
Data *& external_data_pt(const unsigned &i)
Return a pointer to i-th external data object.
void add_internal_data_values_to_vector(Vector< double > &vector_of_values)
Add all internal data and time history values to the vector in the internal storage order.
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
virtual void assign_additional_local_eqn_numbers()
Setup any additional look-up schemes for local equation numbers. Examples of use include using local ...
virtual void get_residuals(Vector< double > &residuals)
Calculate the vector of residuals of the equations in the element. By default initialise the vector t...
virtual void fill_in_contribution_to_dresiduals_dparameter(double *const ¶meter_pt, Vector< double > &dres_dparam)
Add the elemental contribution to the derivatives of the residuals with respect to a parameter....
static bool Suppress_warning_about_any_repeated_data
Static boolean to suppress warnings about repeated data. Defaults to false.
unsigned ninternal_data() const
Return the number of internal data objects.
int Non_halo_proc_ID
Non-halo processor ID for Data; -1 if it's not a halo.
void clear_global_eqn_numbers()
Clear the storage for the global equation numbers and pointers to dofs (if stored)
unsigned Ninternal_data
Number of internal data.
unsigned Nexternal_data
Number of external data.
void add_internal_value_pt_to_map(std::map< unsigned, double * > &map_of_value_pt)
Add pointers to the internal data values to map indexed by the global equation number.
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...
virtual void update_before_internal_fd()
Function that is called before the finite differencing of any internal data. This may be overloaded t...
virtual void fill_in_contribution_to_inner_products(Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
Fill in the contribution to the inner products between given pairs of history values.
virtual void reset_in_internal_fd(const unsigned &i)
Function called within the finite difference loop for internal data after the values in the i-th exte...
virtual void fill_in_contribution_to_djacobian_dparameter(double *const ¶meter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Add the elemental contribution to the derivatives of the elemental Jacobian matrix and residuals with...
virtual void fill_in_contribution_to_hessian_vector_products(Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Fill in contribution to the product of the Hessian (derivative of Jacobian with respect to all variab...
unsigned Ndof
Number of degrees of freedom.
int internal_local_eqn(const unsigned &i, const unsigned &j) const
Return the local equation number corresponding to the j-th value stored at the i-th internal data.
virtual void describe_local_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
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...
virtual void assign_internal_and_external_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers for the internal and external Data This must be called after the gl...
void set_halo(const unsigned &non_halo_proc_ID)
Label the element as halo and specify processor that holds non-halo counterpart.
void add_internal_eqn_numbers_to_vector(Vector< long > &vector_of_eqn_numbers)
Add all equation numbers associated with internal data to the vector in the internal storage order.
virtual void reset_after_internal_fd()
Function that is call after the finite differencing of the internal data. This may be overloaded to r...
virtual unsigned self_test()
Self-test: Have all internal values been classified as pinned/unpinned? Return 0 if OK.
virtual void assign_local_eqn_numbers(const bool &store_local_dof_pt)
Setup the arrays of local equation numbers for the element. If the optional boolean argument is true,...
Data ** Data_pt
Storage for pointers to internal and external data. The data is of the same type and so can be addres...
static bool Suppress_warning_about_repeated_internal_data
Static boolean to suppress warnings about repeated internal data. Defaults to false.
double ** Dof_pt
Storage for array of pointers to degrees of freedom ordered by local equation number....
int external_local_eqn(const unsigned &i, const unsigned &j)
Return the local equation number corresponding to the j-th value stored at the i-th external data.
void describe_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the dofs of the element. The ostream specifies the output stream to which the de...
std::vector< bool > Data_fd
Storage for boolean flags of size Ninternal_data + Nexternal_data that correspond to the data used in...
void flush_external_data()
Flush all external data.
void read_internal_data_values_from_vector(const Vector< double > &vector_of_values, unsigned &index)
Read all internal data and time history values from the vector starting from index....
virtual void update_in_external_fd(const unsigned &i)
Function called within the finite difference loop for external data after a change in any values in t...
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) external data object to the element and return its index (i....
/////////////////////////////////////////////////////////////////////
virtual void dposition_dt(const Vector< double > &zeta, const unsigned &j, Vector< double > &drdt)
j-th time-derivative on object at current time: .
Generic class for numerical integration schemes:
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
double & x(const unsigned &i)
Return the i-th nodal coordinate.
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i). ‘Type’: k; Coordinate direction: i.
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.
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
virtual void node_update(const bool &update_all_time_levels_for_new_node=false)
Interface for functions that update the nodal position using algebraic remeshing strategies....
An OomphLibError object which should be thrown when an run-time error is encountered....
void disable_error_message()
Suppress issueing of the error message in destructor (useful if error is caught successfully!...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
An OomphLibWarning object which should be created as a temporary object to issue a warning....
void shape(const Vector< double > &s, Shape &psi) const
Broken assignment operator.
static PointIntegral Default_integration_scheme
Return spatial dimension of element (=number of local coordinates needed to parametrise the element)
Broken pseudo-integration scheme for points elements: Iit's not clear in general what this integratio...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
unsigned long nindex4() const
Return the range of index 4 of the tensor.
unsigned long nindex2() const
Return the range of index 2 of the tensor.
unsigned long nindex1() const
Return the range of index 1 of the tensor.
unsigned long nindex3() const
Return the range of index 3 of the tensor.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
unsigned lagrangian_dimension() const
Return the number of Lagrangian coordinates that the element requires at all nodes....
void describe_solid_local_dofs(std::ostream &out, const std::string ¤t_string) const
Classifies dofs locally for solid specific aspects.
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.
double raw_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....
int * Position_local_eqn
Array to hold the local equation number information for the solid equations, whatever they may be.
virtual 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...
virtual void assign_solid_local_eqn_numbers(const bool &store_local_dof)
Assigns local equation numbers for the generic solid local equation numbering schemes....
double dshape_lagrangian(const Vector< double > &s, Shape &psi, DShape &dpsidxi) const
Calculate shape functions and derivatives w.r.t. Lagrangian coordinates at local coordinate s....
unsigned nnodal_lagrangian_type() const
Return the number of types of (generalised) nodal Lagrangian coordinates required to interpolate the ...
unsigned Lagrangian_dimension
The Lagrangian dimension of the nodes stored in the element, / i.e. the number of Lagrangian coordina...
virtual void update_before_solid_position_fd()
Function that is called before the finite differencing of any solid position data....
bool Solve_for_consistent_newmark_accel_flag
Flag to indicate which system of equations to solve when assigning initial conditions for time-depend...
SolidInitialCondition * Solid_ic_pt
Pointer to object that specifies the initial condition.
virtual double d2shape_lagrangian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidxi, DShape &d2psidxi) const
Return the geometric shape functions and also first and second derivatives w.r.t. Lagrangian coordina...
void fill_in_generic_jacobian_for_solid_ic(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Helper function to fill in the residuals and (if flag==1) the Jacobian for the setup of an initial co...
virtual 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...
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 double local_to_lagrangian_mapping(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 functi...
virtual 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,...
virtual double interpolated_xi(const Vector< double > &s, const unsigned &i) const
Return i-th FE-interpolated Lagrangian coordinate xi[i] at local coordinate s.
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 interpolated_dxids(const Vector< double > &s, DenseMatrix< double > &dxids) const
Compute derivatives of FE-interpolated Lagrangian coordinates xi with respect to local coordinates: d...
virtual double dshape_lagrangian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidxi) const
Return the geometric shape functions and also first derivatives w.r.t. Lagrangian coordinates at ipt-...
virtual void fill_in_jacobian_from_solid_position_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Use finite differences to calculate the Jacobian entries corresponding to the solid positions....
double multiplier(const Vector< double > &xi)
Access to the "multiplier" for the inertia terms in the consistent determination of the initial condi...
double d2shape_lagrangian(const Vector< double > &s, Shape &psi, DShape &dpsidxi, DShape &d2psidxi) const
Compute the geometric shape functions and also first and second derivatives w.r.t....
void compute_norm(double &el_norm)
Calculate the L2 norm of the displacement u=R-r to overload the compute_norm function in the Generali...
void describe_local_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
virtual ~SolidFiniteElement()
Destructor to clean up any allocated memory.
void fill_in_jacobian_for_newmark_accel(DenseMatrix< double > &jacobian)
Fill in the contributions of the Jacobian matrix for the consistent assignment of the initial "accele...
virtual void reset_after_solid_position_fd()
Function that is call after the finite differencing of the solid position data. This may be overloade...
GeomObject *& geom_object_pt()
(Reference to) pointer to geom object that specifies the initial condition
unsigned & ic_time_deriv()
Which time derivative are we currently assigning?
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
const long & position_eqn_number(const unsigned &k, const unsigned &i) const
Return the equation number for generalised Eulerian coordinate: type of coordinate: k,...
void describe_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the dofs of the Node. The ostream specifies the output stream to which the descr...
Data *const & variable_position_pt() const
Pointer to variable_position data (const version)
////////////////////////////////////////////////////////////////////// //////////////////////////////...
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
std::string type() const
Return string that indicates the type of the timestepper (e.g. "BDF", "Newmark", etc....
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned Max_newton_iterations
Maximum number of newton iterations.
double Newton_tolerance
Convergence tolerance for the newton solver.
unsigned N_local_points
Number of points along one dimension of each element used to create coordinates within the element in...
double Radius_multiplier_for_fast_exit_from_locate_zeta
Multiplier for (zeta-based) outer radius of element used for deciding that point is outside element....
const double Pi
50 digits from maple
double dot(const Vector< double > &a, const Vector< double > &b)
Probably not always best/fastest because not optimised for dimension but useful...
void cross(const Vector< double > &A, const Vector< double > &B, Vector< double > &C)
Cross product using "proper" output (move semantics means this is ok nowadays).
double angle(const Vector< double > &a, const Vector< double > &b)
Get the angle between two vector.
double magnitude(const Vector< double > &a)
Get the magnitude of a vector.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...