Commit 2591d8be authored by Guido Giuntoli's avatar Guido Giuntoli

Adding geometrical function for fibers.

Deleting useless function mvp_3 in util.h.
Overloading mvp and replacing mvp_2 & mvp_3.
parent 2e46d29f
......@@ -72,7 +72,7 @@ inline void print_vec(const double *vec, int n, const char file_name[])
}
inline void mvp_2(const double m[2][2], const double x[2], double y[2])
inline void mvp(const double m[2][2], const double x[2], double *y)
{
for (int i = 0; i < 2; ++i) {
double tmp = 0.0;
......@@ -83,25 +83,67 @@ inline void mvp_2(const double m[2][2], const double x[2], double y[2])
}
inline void mvp_3(const double m[2][2], const double x[2], double y[2])
inline void mvp(const double m[3][3], const double x[3], double *y)
{
for (int i = 0; i < 2; ++i) {
for (int i = 0; i < 3; ++i) {
double tmp = 0.0;
for (int j = 0; j < 2; ++j)
for (int j = 0; j < 3; ++j)
tmp += m[i][j] * x[j];
y[i] = tmp;
}
}
inline void mvp_3(const double m[3][3], const double x[3], double y[3])
template<typename T, int n>
inline double norm(const T v1[n])
{
for (int i = 0; i < 3; ++i) {
double tmp = 0.0;
for (int j = 0; j < 3; ++j)
tmp += m[i][j] * x[j];
y[i] = tmp;
/*
* Returns sqrt (v1[0] * v1[0] + ... + v1[n-1] * v1[n-1])
*/
T tmp = 0;
for (int i = 0; i < n; ++i) {
tmp += v1[i] * v1[i];
}
return sqrt((double)tmp);
}
template<typename T, int n>
inline T dot_prod(const T v1[n], const T v2[n])
{
/*
* Returns v1[0] * v2[0] + ... + v1[n-1] * v2[n-1]
*/
T tmp = 0;
for (int i = 0; i < n; ++i) {
tmp += v1[i] * v2[i];
}
return tmp;
}
inline bool point_inf_cilinder(const double dir[3], const double center[3],
const double radius, const double point[3])
{
/*
* Returns <true> if <point> is inside the infinite cilinder with
* <direction>, <center> and <radius>. Returns <false> otherwise.
*/
const double v[3] = {
point[0] - center[0],
point[1] - center[1],
point[2] - center[2]};
const double dir_dot_v = dot_prod<double, 3>(dir, v);
const double norm_dir = norm<double, 3>(dir);
const double norm_v = norm<double, 3>(v);
const double cos_tetha = dir_dot_v / (norm_dir * norm_v);
const double sin_tetha = sqrt(1 - cos_tetha * cos_tetha);
const double d_1 = norm_v * sin_tetha;
return (d_1 <= radius) ? true : false;
}
......
......@@ -35,25 +35,25 @@ void micropp<2>::set_displ_bc(const double eps[nvoi], double *u)
for (int i = 0; i < nx; ++i) {
const int n = nod_index(i, 0, 0); // y = 0
const double coor[2] = { i * dx, 0 };
mvp_2(eps_t, coor, &u[n * dim]);
mvp(eps_t, coor, &u[n * dim]);
}
for (int i = 0; i < nx; ++i) {
const int n = nod_index(i, ny - 1, 0); // y = ly
const double coor[2] = { i * dx, ly };
mvp_2(eps_t, coor, &u[n * dim]);
mvp(eps_t, coor, &u[n * dim]);
}
for (int j = 1; j < ny - 1; ++j) {
const int n = nod_index(0, j, 0); // x = 0
const double coor[2] = { 0, j * dy };
mvp_2(eps_t, coor, &u[n * dim]);
mvp(eps_t, coor, &u[n * dim]);
}
for (int j = 1; j < ny - 1; ++j) {
const int n = nod_index(nx - 1, j, 0); // x = lx
const double coor[2] = { lx, j * dy };
mvp_2(eps_t, coor, &u[n * dim]);
mvp(eps_t, coor, &u[n * dim]);
}
}
......
......@@ -37,7 +37,7 @@ void micropp<3>::set_displ_bc(const double eps[nvoi], double *u)
for (int j = 0; j < ny; ++j) {
const int n = nod_index3D(i, j, 0); // z = 0
const double coor[3] = { i * dx, j * dy, 0 };
mvp_3(eps_t, coor, &u[n * dim]);
mvp(eps_t, coor, &u[n * dim]);
}
}
......@@ -45,7 +45,7 @@ void micropp<3>::set_displ_bc(const double eps[nvoi], double *u)
for (int j = 0; j < ny; ++j) {
const int n = nod_index3D(i, j, nz - 1); // z = lz
const double coor[3] = { i * dx, j * dy, lz };
mvp_3(eps_t, coor, &u[n * dim]);
mvp(eps_t, coor, &u[n * dim]);
}
}
......@@ -53,7 +53,7 @@ void micropp<3>::set_displ_bc(const double eps[nvoi], double *u)
for (int k = 0; k < nz; ++k) {
const int n = nod_index3D(i, 0, k); // y = 0
const double coor[3] = { i * dx, 0, k * dz };
mvp_3(eps_t, coor, &u[n * dim]);
mvp(eps_t, coor, &u[n * dim]);
}
}
......@@ -61,7 +61,7 @@ void micropp<3>::set_displ_bc(const double eps[nvoi], double *u)
for (int k = 0; k < nz; ++k) {
const int n = nod_index3D(i, ny - 1, k); // y = ly
const double coor[3] = { i * dx, ly , k * dz };
mvp_3(eps_t, coor, &u[n * dim]);
mvp(eps_t, coor, &u[n * dim]);
}
}
......@@ -69,7 +69,7 @@ void micropp<3>::set_displ_bc(const double eps[nvoi], double *u)
for (int k = 0; k < nz; ++k) {
const int n = nod_index3D(0, j, k); // x = 0
const double coor[3] = { 0, j * dy , k * dz };
mvp_3(eps_t, coor, &u[n * dim]);
mvp(eps_t, coor, &u[n * dim]);
}
}
......@@ -77,7 +77,7 @@ void micropp<3>::set_displ_bc(const double eps[nvoi], double *u)
for (int k = 0; k < nz; ++k) {
const int n = nod_index3D(nx - 1, j, k); // x = lx
const double coor[3] = { lx, j * dy , k * dz };
mvp_3(eps_t, coor, &u[n * dim]);
mvp(eps_t, coor, &u[n * dim]);
}
}
}
......
......@@ -39,7 +39,7 @@ int main (int argc, char *argv[])
const double s2_exact[2] = { 532.3834, -246.7034 };
double s2[2];
mvp_2(a1, v1, s2);
mvp(a1, v1, s2);
for (int i = 0; i < 2; ++i)
assert(fabs(s2[i] - s2_exact[i]) < 1.0e-10);
......@@ -51,7 +51,7 @@ int main (int argc, char *argv[])
const double s3_exact[3] = { -126.4282, -53.53220, -154.4488 };
double s3[3];
mvp_3(a3, v3, s3);
mvp(a3, v3, s3);
for (int i = 0; i < 3; ++i)
assert(fabs(s3[i] - s3_exact[i]) < 1.0e-10);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment