+inline void Math::matrixEig(__attribute__((unused)) UInt n,
+ __attribute__((unused)) T * A,
+ __attribute__((unused)) T * d,
+ __attribute__((unused)) T * V) {
+ AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function! Or implement it for the type " << debug::demangle(typeid(T).name()));
+}
+
#ifdef AKANTU_USE_LAPACK
-inline void Math::matrixEig(UInt n, Real * A, Real * d, Real * V) {
+template<>
+inline void Math::matrixEig<double>(UInt n, double * A, double * d, double * V) {
// Matrix A is row major, so the lapack function in fortran will process
// A^t. Asking for the left eigenvectors of A^t will give the transposed right
// eigenvectors of A so in the C++ code the right eigenvectors.
char jobvl;
if(V != NULL)
jobvl = 'V'; // compute left eigenvectors
else
jobvl = 'N'; // compute left eigenvectors
char jobvr('N'); // compute right eigenvectors
double * di = new double[n]; // imaginary part of the eigenvalues
+inline void Math::inv(__attribute__((unused)) UInt n,
+ __attribute__((unused)) const T * A,
+ __attribute__((unused)) T * Ainv) {
+ AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function! Or implement it for the type " << debug::demangle(typeid(T).name()));
+}
+
+#ifdef AKANTU_USE_LAPACK
+template<>
+inline void Math::inv<double>(UInt n, const double * A, double * invA) {
+ int N = n;
+ int info;
+ int * ipiv = new int[N+1];
+ int lwork = N*N;
+ double * work = new double[lwork];
+
+ std::copy(A, A + n*n, invA);
+
+ dgetrf_(&N, &N, invA, &N, ipiv, &info);
+ if(info > 0) {
+ AKANTU_DEBUG_ERROR("Singular matrix - cannot factorize it (info: "
+ * interpolate the field on a point given in natural coordinates the field which
+ * values are given on the node of the element
+ *
+ * @param natural_coords natural coordinates of point where to interpolate \xi
+ * @param nodal_values values of the function per node @f$ f_{ij} = f_{n_i j} @f$ so it should be a matrix of size nb_nodes_per_element @f$\times@f$ nb_degree_of_freedom
+ * @param interpolated interpolated value of f @f$ f_j(\xi) = \sum_i f_{n_i j} N_i @f$
std::cout << " * Historic Velocity Friction Coefficient: total history time = " << total_time << " with numbers of time steps = " << nb_time_steps << std::endl;
// declare vector
this->weights = new Vector<Real>(nb_time_steps,1);
this->historic_velocities = new CircularVector< Vector<Real> * >(nb_time_steps, 1);
// fill the weights vector
Real time = 0.;
for (Int i=nb_time_steps-1; i>=0; --i) {
(*weights)(i) = exp(-beta*time);
time += time_step;
}
this->generalized_sliding_velocities = new Vector<Real>(0,1);
// find index of master surface in impactors_information
Int master_index = -1;
for (UInt i=0; i < impactor_information_vector.size(); ++i) {
if (impactor_information_vector.at(i)->master_id == this->master_surface) {
master_index = i;
break;
}
}
AKANTU_DEBUG_ASSERT(master_index != -1, "No impactor information object for master" << this->master_surface << "in impactors_information vector that is ");
AKANTU_DEBUG_INFO("Receiving list of elements (" << status.getSize() - 1 << " elements) no longer needed by proc " << p << " TAG("<< Tag::genTag(p, 0, 0) <<")");
comm.receive(list.storage(), list.getSize(),
p, Tag::genTag(p, 0, 0));
if(list.getSize() == 1 && list(0) == 0) continue;
list.erase(list.getSize() - 1);
Vector<Element> new_send;
for (UInt ns = 0; ns < list.getSize(); ++ns) {
new_send.push_back(send(list(ns)));
}
AKANTU_DEBUG_INFO("I had " << send.getSize() << " elements to send to proc " << p << " and "