commit 095d8ea2ae62962f45a6557b19e0abf542fbead9
Author: Julio C. Ramirez-Ceballos <julio@jcrcx.xyz>
Date: Wed, 27 Nov 2024 11:37:34 +0100
Initial commit
Diffstat:
| A | Makefile | | | 2 | ++ |
| A | clt.h | | | 788 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | laminate.dat | | | 7 | +++++++ |
| A | main | | | 0 | |
| A | main.c | | | 157 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | materials.dat | | | 14 | ++++++++++++++ |
| A | rw.c | | | 218 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | utils.h | | | 260 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
8 files changed, 1446 insertions(+), 0 deletions(-)
diff --git a/Makefile b/Makefile
@@ -0,0 +1,2 @@
+all:
+ gcc -Wall main.c -lm -o main
diff --git a/clt.h b/clt.h
@@ -0,0 +1,788 @@
+/*
+ * CLASSICAL LAMINATE THEORY
+ * clt.h
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "utils.h"
+#include "rw.c"
+
+#define MAX_LAYERS 100 /* Assume a maximum number of layers for simplicity */
+#define MAX_MODE_LENGTH 10
+#define NELEMS(x) (sizeof(x) / sizeof((x)[0]))
+#define N 6 // Example for a 6x6 matrix, adjust for your needs
+
+/* Custom error for laminate layup issues */
+typedef struct {
+ char *message;
+} LaminateLayupError;
+
+typedef struct {
+ double *stress_inf, *stress_sup, *strain_inf, *strain_sup;
+} StressStrainResult;
+
+typedef struct {
+ double fs;
+ char *mode;
+} FailureResult;
+
+typedef struct {
+ char mode[MAX_MODE_LENGTH];
+ double fs;
+} FSResult;
+
+/*void calc_stressCLT(MaterialProperties *mat_list, Laminate *lam, */
+/* double F[6], int *fail_list, double dT, double dM);*/
+void assemble_Z(double *Z, Laminate *lam);
+void calc_thermal_forces(double *Nt, MaterialProperties *mat_list,
+ Laminate *lam, double *Z, int *fail_list, double dT);
+void calc_moisture_forces(double *Nm, MaterialProperties *mat_list,
+ Laminate *lam, double *Z, int *fail_list, double dM);
+void assemble_Q(double Q[3][3], MaterialProperties *mat_prop,
+ int fail_type);
+void assemble_T(double T[3][3], double angle);
+void assemble_ABD(double ABD[6][6], MaterialProperties *mat_list,
+ Laminate *lam, double *Z, int *fail_list);
+void fs_hashin_2D(MaterialProperties *mat_prop, double sig1, double sig2,
+ double tau, FSResult *result);
+void hashin_2D(MaterialProperties *mat_list, Laminate *lam,
+ double **stress_inf, double **stress_sup,
+ FSResult *fs_inf, FSResult *fs_sup);
+void fs_tsaiwu_2D(MaterialProperties *mat_prop,
+ double sig1, double sig2, double tau,
+ FSResult *result);
+void tsaiwu_2D(MaterialProperties *mat_list, Laminate *lam,
+ double **stress_inf, double **stress_sup,
+ FSResult *fs_inf, FSResult *fs_sup);
+void fs_maxstress_2D(MaterialProperties *mat_prop, double sig1, double sig2,
+ double tau, FSResult *result);
+void maxstress_2D(MaterialProperties *mat_list, Laminate *lam,
+ double **stress_inf, double **stress_sup,
+ FSResult *fs_inf, FSResult *fs_sup);
+void print_material_properties(MaterialProperties *mat, int mat_id);
+
+void
+assemble_Z(double *Z, Laminate *lam)
+{
+ int i;
+ double total_thk = 0.0;
+
+ for (i = 0; i < lam->num_layers; ++i) {
+ total_thk += lam->thk[i];
+ }
+
+ Z[0] = -total_thk / 2.0;
+ for (i = 1; i <= lam->num_layers; ++i) {
+ Z[i] = Z[i-1] + lam->thk[i-1];
+ }
+}
+
+void
+calc_thermal_forces(double *Nt, MaterialProperties *mat_list,
+ Laminate *lam, double *Z, int *fail_list, double dT)
+{
+ double a[3], a_LCS[3];
+ double Q[3][3], Qbar[3][3], T[3][3], Ti[3][3], temp[3][3];
+ int i, j, k;
+
+ for (i = 0; i < lam->num_layers; ++i) {
+ MaterialProperties *mat_prop = &mat_list[lam->mat_id[i]];
+ a[0] = mat_prop->a1; a[1] = mat_prop->a2; a[2] = 0.0;
+
+ assemble_T(T, lam->ang[i]);
+
+ /* Transform to Laminate Coordinate System */
+ for (j = 0; j < 3; ++j) {
+ for (k = 0; k < 3; ++k) {
+ a_LCS[j] += T[k][j] * a[k];
+ }
+ }
+
+ /* Assemble Q matrix */
+ assemble_Q(Q, mat_prop, fail_list ? fail_list[i] : 0);
+
+ /* Calculate Qbar */
+ transpose(T,Ti);
+ matmul(Ti,Q,temp);
+ matmul(Q,T,Qbar);
+
+ /* Calculate force resultants */
+ for (j = 0; j < 3; ++j) {
+ for (k = 0; k < 3; ++k) {
+ double force = Qbar[j][k] * a_LCS[k] * dT;
+ Nt[j] += force * (Z[i+1] - Z[i]);
+ Nt[j+3] += force * 0.5 *
+ (Z[i+1]*Z[i+1] - Z[i]*Z[i]);
+ }
+ }
+ }
+}
+
+void
+calc_moisture_forces(double *Nm, MaterialProperties *mat_list,
+ Laminate *lam, double *Z, int *fail_list, double dM)
+{
+ double b[3], b_LCS[3];
+ double Q[3][3], Qbar[3][3], T[3][3], Ti[3][3], temp[3][3];
+ int i, j, k;
+
+ for (i = 0; i < lam->num_layers; ++i) {
+ MaterialProperties *mat_prop = &mat_list[lam->mat_id[i]];
+ b[0] = mat_prop->b1; b[1] = mat_prop->b2; b[2] = 0.0;
+
+ /* Apply transformed matrix */
+ assemble_T(T, lam->ang[i]);
+
+ /* Transform to Laminate Coordinate System */
+ for (j = 0; j < 3; ++j) {
+ for (k = 0; k < 3; ++k) {
+ b_LCS[j] += T[k][j] * b[k];
+ }
+ }
+
+ /* Assemble Q matrix */
+ assemble_Q(Q, mat_prop, fail_list ? fail_list[i] : 0);
+
+ /* Calculate Qbar */
+ transpose(T,Ti);
+ matmul(Ti,Q,temp);
+ matmul(Q,T,Qbar);
+
+ /* Calculate force resultants */
+ for (j = 0; j < 3; ++j) {
+ for (k = 0; k < 3; ++k) {
+ double force = Qbar[j][k] * b_LCS[k] * dM;
+ Nm[j] += force * (Z[i+1] - Z[i]);
+ Nm[j+3] += force * 0.5 *
+ (Z[i+1]*Z[i+1] - Z[i]*Z[i]);
+ }
+ }
+ }
+}
+
+void
+assemble_Q(double Q[3][3], MaterialProperties *mat_prop, int fail_type)
+{
+
+ if (mat_prop->isotropic) {
+ double E, nu, E_div;
+ E = mat_prop->E1;
+ nu = mat_prop->nu12;
+
+ if (fail_type == 1 || fail_type == 2 || fail_type == 3) {
+ double df = 0.001;
+ E *= df;
+ nu *= df;
+ }
+
+ E_div = E / (1 - nu * nu);
+
+ Q[0][0] = Q[1][1] = E_div;
+ Q[0][1] = Q[1][0] = nu * E_div;
+ Q[2][2] = E_div * (1 - nu) / 2.0; // Since G = E / (2(1+nu)) for isotropic
+ Q[0][2] = Q[1][2] = Q[2][0] = Q[2][1] = 0.0;
+ } else {
+ double E1 = mat_prop->E1,
+ E2 = mat_prop->E2,
+ nu12 = mat_prop->nu12,
+ G12 = mat_prop->G12;
+ double df = 0.001; /* Degradation factor */
+
+ if (fail_type == 1 || fail_type == 2) { /* fiber or shear failure */
+ E1 *= df;
+ E2 *= df;
+ nu12 *= df;
+ G12 *= df;
+ } else if (fail_type == 3) { /* matrix failure */
+ E2 *= df;
+ nu12 *= df;
+ G12 *= df;
+ }
+
+ double n21 = nu12 * E2 / E1;
+
+ Q[0][0] = E1 / (1 - nu12 * n21);
+ Q[0][1] = Q[1][0] = nu12 * E1 * E2 / (E1 - nu12 * nu12 * E2);
+ Q[1][1] = E1 * E2 / (E1 - nu12 * nu12 * E2);
+ Q[2][2] = G12;
+ Q[0][2] = Q[1][2] = Q[2][0] = Q[2][1] = 0.0;
+ }
+}
+
+void
+assemble_T(double T[3][3], double angle)
+{
+ double rad = angle * DEGTORAD;
+
+ double cos = Cos(rad), sin = Sin(rad);
+ double cs = cos * sin, cc = cos * cos, ss = sin * sin;
+
+ T[0][0] = cc; T[0][1] = ss; T[0][2] = cs;
+ T[1][0] = ss; T[1][1] = cc; T[1][2] = -cs;
+ T[2][0] = -2*cs; T[2][1] = 2*cs; T[2][2] = cc - ss;
+}
+
+void
+assemble_ABD(double ABD[6][6], MaterialProperties *mat_list,
+ Laminate *lam, double *Z, int *fail_list)
+{
+ double A[3][3] = {0}, B[3][3] = {0}, D[3][3] = {0};
+ double Q[3][3], T[3][3], Ti[3][3];
+ double temp[3][3], Qbar[3][3];
+
+ for (int i = 0; i < lam->num_layers; ++i) {
+ MaterialProperties *mat_prop = &mat_list[lam->mat_id[i]];
+
+ assemble_Q(Q, mat_prop, fail_list ? fail_list[i] : 0);
+ assemble_T(T, lam->ang[i]);
+
+ /* Transpose T to get Ti */
+ transpose(T,Ti);
+
+ /* Qbar = T^T * Q * T */
+ matmul(Ti, Q, temp);
+ matmul(temp, T, Qbar);
+
+ /* B and D matrices */
+ //double h = Z[i+1] - Z[i]; /* thickness of the layer */
+ for (int j = 0; j < 3; ++j) {
+ for (int k = 0; k < 3; ++k) {
+ A[j][k] += Qbar[j][k] * (Z[i+1] - Z[i]);
+ B[j][k] += 0.5 * Qbar[j][k] * (Z[i+1] * Z[i+1] - Z[i] * Z[i]);
+ D[j][k] += CONST13 * Qbar[j][k] * (Z[i+1] * Z[i+1] * Z[i+1] - Z[i] * Z[i] * Z[i]);
+ }
+ }
+ }
+
+ /* Combine into ABD matrix */
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ ABD[i][j] = A[i][j]; /* A matrix */
+ ABD[i][j+3] = ABD[i+3][j] = B[i][j]; /* B matrix */
+ ABD[i+3][j+3] = D[i][j]; /* D matrix */
+ }
+ }
+}
+
+/* Tsai-Wu criterion */
+void
+fs_tsaiwu_2D(MaterialProperties *mat_prop,
+ double sig1, double sig2, double tau,
+ FSResult *result)
+{
+ double f11 = 1 / (mat_prop->Xt * mat_prop->Xc);
+ double f22 = 1 / (mat_prop->Yt * mat_prop->Yc);
+ double f12 = -1 / (2 * sqrt(mat_prop->Xt * mat_prop->Xc * mat_prop->Yt * mat_prop->Yc));
+ double f66 = 1 / (mat_prop->S12 * mat_prop->S12);
+ double f1 = 1/mat_prop->Xt - 1/mat_prop->Xc;
+ double f2 = 1/mat_prop->Yt - 1/mat_prop->Yc;
+
+ double a = f11*sig1*sig1 + f22*sig2*sig2 + f66*tau*tau + 2*f12*sig1*sig2;
+ double b = f1*sig1 + f2*sig2;
+
+ result->fs = (-b + sqrt(b*b + 4*a)) / (2*a);
+
+ double H1 = fabs(f1*sig1 + f11*sig1*sig1);
+ double H2 = fabs(f2*sig2 + f22*sig2*sig2);
+ double H6 = fabs(f66*tau*tau);
+
+ if (H1 >= H2 && H1 >= H6) {
+ strcpy(result->mode, "fiber");
+ } else if (H2 >= H1 && H2 >= H6) {
+ strcpy(result->mode, "matrix");
+ } else {
+ strcpy(result->mode, "shear");
+ }
+}
+
+void
+tsaiwu_2D(MaterialProperties *mat_list, Laminate *lam, double **stress_inf,
+ double **stress_sup, FSResult *fs_inf, FSResult *fs_sup)
+{
+ for (int i = 0; i < lam->num_layers; i++) {
+ MaterialProperties *mat_prop = &mat_list[lam->mat_id[i]];
+
+ /* Superior Stresses */
+ fs_tsaiwu_2D(mat_prop, stress_sup[0][i], stress_sup[1][i], stress_sup[2][i], &fs_sup[i]);
+
+ /* Inferior Stresses */
+ fs_tsaiwu_2D(mat_prop, stress_inf[0][i], stress_inf[1][i], stress_inf[2][i], &fs_inf[i]);
+ }
+}
+
+/* Maximum Stress criterion */
+void
+fs_maxstress_2D(MaterialProperties *mat_prop, double sig1, double sig2,
+ double tau, FSResult *result)
+{
+ double f_1 = sig1 > 0 ? sig1 / mat_prop->Xt : -sig1 / mat_prop->Xc;
+ double f_2 = sig2 > 0 ? sig2 / mat_prop->Yt : -sig2 / mat_prop->Yc;
+ double f_s = fabs(tau) / mat_prop->S12;
+
+ double f_max = fmax(fmax(f_1, f_2), f_s);
+
+ result->fs = 1 / f_max;
+
+ if (f_max == f_1) {
+ strcpy(result->mode, "fiber");
+ } else if (f_max == f_2) {
+ strcpy(result->mode, "matrix");
+ } else {
+ strcpy(result->mode, "shear");
+ }
+}
+
+void
+maxstress_2D(MaterialProperties *mat_list, Laminate *lam,
+ double **stress_inf, double **stress_sup,
+ FSResult *fs_inf, FSResult *fs_sup)
+{
+ for (int i = 0; i < lam->num_layers; i++) {
+ MaterialProperties *mat_prop = &mat_list[lam->mat_id[i]];
+
+ /* Superior Stresses */
+ fs_maxstress_2D(mat_prop, stress_sup[0][i], stress_sup[1][i],
+ stress_sup[2][i], &fs_sup[i]);
+
+ /* Inferior Stresses */
+ fs_maxstress_2D(mat_prop, stress_inf[0][i], stress_inf[1][i],
+ stress_inf[2][i], &fs_inf[i]);
+ }
+}
+
+/* Maximum Strain criterion is not directly translated due to lack of strain input in the given code. */
+
+/* Hashin criterion */
+void
+fs_hashin_2D(MaterialProperties *mat_prop,
+ double sig1, double sig2, double tau,
+ FSResult *result)
+{
+ double f_1 = sig1 >= 0 ? sig1 / mat_prop->Xt : -sig1 / mat_prop->Xc;
+ double f_2 = sig2 >= 0 ? sqrt((sig2/mat_prop->Yt)*(sig2/mat_prop->Yt)
+ + (tau/mat_prop->S12)*(tau/mat_prop->S12)) :
+ sqrt((sig2/mat_prop->Yc)*(sig2/mat_prop->Yc) +
+ (tau/mat_prop->S12)*(tau/mat_prop->S12));
+
+ double f_max = fmax(f_1, f_2);
+
+ result->fs = 1 / f_max;
+
+ if (f_max == f_1) {
+ strcpy(result->mode, "fiber");
+ } else {
+ strcpy(result->mode, "matrix");
+ }
+}
+
+void
+hashin_2D(MaterialProperties *mat_list, Laminate *lam,
+ double **stress_inf, double **stress_sup,
+ FSResult *fs_inf, FSResult *fs_sup)
+{
+ for (int i = 0; i < lam->num_layers; i++) {
+ MaterialProperties *mat_prop = &mat_list[lam->mat_id[i]];
+
+ /* Superior Stresses */
+ fs_hashin_2D(mat_prop, stress_sup[0][i], stress_sup[1][i], stress_sup[2][i], &fs_sup[i]);
+
+ /* Inferior Stresses */
+ fs_hashin_2D(mat_prop, stress_inf[0][i], stress_inf[1][i], stress_inf[2][i], &fs_inf[i]);
+ }
+}
+
+double E_x(double ABD[6][6], double h) {
+ double matrix1[6][6], matrix2[5][5];
+ double det1, det2, Ex;
+
+ // Copy relevant parts of ABD to matrix1 for the first determinant
+ for (int i = 0; i < 6; i++) {
+ for (int j = 0; j < 6; j++) {
+ matrix1[i][j] = ABD[i][j];
+ }
+ }
+
+ // Perform LU decomposition for matrix1
+ if (LUDecompose(&matrix1[0][0], 6, 6) < 0) {
+ fprintf(stderr, "Error during LU decomposition for matrix1.\n");
+ return -1;
+ }
+
+ // Calculate determinant for matrix1
+ det1 = LUDeterminant(&matrix1[0][0], 6, 6);
+
+ matrix2[0][0] = ABD[1][1];
+ matrix2[0][1] = ABD[1][2];
+ matrix2[0][2] = ABD[1][3];
+ matrix2[0][3] = ABD[1][4];
+ matrix2[0][4] = ABD[1][5];
+ matrix2[1][0] = ABD[2][1];
+ matrix2[1][1] = ABD[2][2];
+ matrix2[1][2] = ABD[2][3];
+ matrix2[1][3] = ABD[2][4];
+ matrix2[1][4] = ABD[2][5];
+ matrix2[2][0] = ABD[3][1];
+ matrix2[2][1] = ABD[3][2];
+ matrix2[2][2] = ABD[3][3];
+ matrix2[2][3] = ABD[3][4];
+ matrix2[2][4] = ABD[3][5];
+ matrix2[3][0] = ABD[4][1];
+ matrix2[3][1] = ABD[4][2];
+ matrix2[3][2] = ABD[4][3];
+ matrix2[3][3] = ABD[4][4];
+ matrix2[3][4] = ABD[4][5];
+ matrix2[4][0] = ABD[5][1];
+ matrix2[4][1] = ABD[5][2];
+ matrix2[4][2] = ABD[5][3];
+ matrix2[4][3] = ABD[5][4];
+ matrix2[4][4] = ABD[5][5];
+
+ // Perform LU decomposition for matrix2
+ if (LUDecompose(&matrix2[0][0], 5, 5) < 0) {
+ fprintf(stderr, "Error during LU decomposition for matrix2.\n");
+ return -1;
+ }
+
+ // Calculate determinant for matrix2
+ det2 = LUDeterminant(&matrix2[0][0], 5, 5);
+
+ // Check for division by zero
+ if (det2 == 0) {
+ fprintf(stderr, "Error: Division by zero in E_x calculation, determinant of the submatrix is zero.\n");
+ return -1;
+ }
+
+ // Compute E_x
+ Ex = (det1 / det2) / h;
+ return Ex;
+}
+
+double E_y(double ABD[6][6], double h) {
+ double matrix1[6][6], matrix2[5][5];
+ double det1, det2, Ey;
+
+ // Copy relevant parts of ABD to matrix1 for the first determinant
+ for (int i = 0; i < 6; i++) {
+ for (int j = 0; j < 6; j++) {
+ matrix1[i][j] = ABD[i][j];
+ }
+ }
+
+ // Perform LU decomposition for matrix1
+ if (LUDecompose(&matrix1[0][0], 6, 6) < 0) {
+ fprintf(stderr, "Error during LU decomposition for matrix1.\n");
+ return -1;
+ }
+
+ // Calculate determinant for matrix1
+ det1 = LUDeterminant(&matrix1[0][0], 6, 6);
+
+ matrix2[0][0] = ABD[0][0];
+ matrix2[0][1] = ABD[0][2];
+ matrix2[0][2] = ABD[0][3];
+ matrix2[0][3] = ABD[0][4];
+ matrix2[0][4] = ABD[0][5];
+ matrix2[1][0] = ABD[2][0];
+ matrix2[1][1] = ABD[2][2];
+ matrix2[1][2] = ABD[2][3];
+ matrix2[1][3] = ABD[2][4];
+ matrix2[1][4] = ABD[2][5];
+ matrix2[2][0] = ABD[3][0];
+ matrix2[2][1] = ABD[3][2];
+ matrix2[2][2] = ABD[3][3];
+ matrix2[2][3] = ABD[3][4];
+ matrix2[2][4] = ABD[3][5];
+ matrix2[3][0] = ABD[4][0];
+ matrix2[3][1] = ABD[4][2];
+ matrix2[3][2] = ABD[4][3];
+ matrix2[3][3] = ABD[4][4];
+ matrix2[3][4] = ABD[4][5];
+ matrix2[4][0] = ABD[5][0];
+ matrix2[4][1] = ABD[5][2];
+ matrix2[4][2] = ABD[5][3];
+ matrix2[4][3] = ABD[5][4];
+ matrix2[4][4] = ABD[5][5];
+
+ // Perform LU decomposition for matrix2
+ if (LUDecompose(&matrix2[0][0], 5, 5) < 0) {
+ fprintf(stderr, "Error during LU decomposition for matrix2.\n");
+ return -1;
+ }
+
+ // Calculate determinant for matrix2
+ det2 = LUDeterminant(&matrix2[0][0], 5, 5);
+
+ // Check for division by zero
+ if (det2 == 0) {
+ fprintf(stderr, "Error: Division by zero in E_x calculation, determinant of the submatrix is zero.\n");
+ return -1;
+ }
+
+ // Compute E_x
+ Ey = (det1 / det2) / h;
+ return Ey;
+}
+
+double G_xy(double ABD[6][6], double h) {
+ double matrix1[6][6], matrix2[5][5];
+ double det1, det2, Gxy;
+
+ // Copy relevant parts of ABD to matrix1 for the first determinant
+ for (int i = 0; i < 6; i++) {
+ for (int j = 0; j < 6; j++) {
+ matrix1[i][j] = ABD[i][j];
+ }
+ }
+
+ // Perform LU decomposition for matrix1
+ if (LUDecompose(&matrix1[0][0], 6, 6) < 0) {
+ fprintf(stderr, "Error during LU decomposition for matrix1.\n");
+ return -1;
+ }
+
+ // Calculate determinant for matrix1
+ det1 = LUDeterminant(&matrix1[0][0], 6, 6);
+
+ matrix2[0][0] = ABD[0][0];
+ matrix2[0][1] = ABD[0][1];
+ matrix2[0][2] = ABD[0][3];
+ matrix2[0][3] = ABD[0][4];
+ matrix2[0][4] = ABD[0][5];
+
+ matrix2[1][0] = ABD[1][0];
+ matrix2[1][1] = ABD[1][1];
+ matrix2[1][2] = ABD[1][3];
+ matrix2[1][3] = ABD[1][4];
+ matrix2[1][4] = ABD[1][5];
+
+ matrix2[2][0] = ABD[3][0];
+ matrix2[2][1] = ABD[3][1];
+ matrix2[2][2] = ABD[3][3];
+ matrix2[2][3] = ABD[3][4];
+ matrix2[2][4] = ABD[3][5];
+
+ matrix2[3][0] = ABD[4][0];
+ matrix2[3][1] = ABD[4][1];
+ matrix2[3][2] = ABD[4][3];
+ matrix2[3][3] = ABD[4][4];
+ matrix2[3][4] = ABD[4][5];
+
+ matrix2[4][0] = ABD[5][0];
+ matrix2[4][1] = ABD[5][1];
+ matrix2[4][2] = ABD[5][3];
+ matrix2[4][3] = ABD[5][4];
+ matrix2[4][4] = ABD[5][5];
+
+ // Perform LU decomposition for matrix2
+ if (LUDecompose(&matrix2[0][0], 5, 5) < 0) {
+ fprintf(stderr, "Error during LU decomposition for matrix2.\n");
+ return -1;
+ }
+
+ // Calculate determinant for matrix2
+ det2 = LUDeterminant(&matrix2[0][0], 5, 5);
+
+ // Check for division by zero
+ if (det2 == 0) {
+ fprintf(stderr, "Error: Division by zero in E_x calculation, determinant of the submatrix is zero.\n");
+ return -1;
+ }
+
+ // Compute E_x
+ Gxy = (det1 / det2) / h;
+ return Gxy;
+}
+
+double nu_xy(double ABD[6][6], double h) {
+ double matrix1[5][5], matrix2[5][5];
+ double det1, det2, nxy;
+
+ matrix1[0][0] = ABD[1][0];
+ matrix1[0][1] = ABD[1][2];
+ matrix1[0][2] = ABD[1][3];
+ matrix1[0][3] = ABD[1][4];
+ matrix1[0][4] = ABD[1][5];
+
+ matrix1[1][0] = ABD[2][0];
+ matrix1[1][1] = ABD[2][2];
+ matrix1[1][2] = ABD[2][3];
+ matrix1[1][3] = ABD[2][4];
+ matrix1[1][4] = ABD[2][5];
+
+ matrix1[2][0] = ABD[3][0];
+ matrix1[2][1] = ABD[3][2];
+ matrix1[2][2] = ABD[3][3];
+ matrix1[2][3] = ABD[3][4];
+ matrix1[2][4] = ABD[3][5];
+
+ matrix1[3][0] = ABD[4][0];
+ matrix1[3][1] = ABD[4][2];
+ matrix1[3][2] = ABD[4][3];
+ matrix1[3][3] = ABD[4][4];
+ matrix1[3][4] = ABD[4][5];
+
+ matrix1[4][0] = ABD[5][0];
+ matrix1[4][1] = ABD[5][2];
+ matrix1[4][2] = ABD[5][3];
+ matrix1[4][3] = ABD[5][4];
+ matrix1[4][4] = ABD[5][5];
+
+ /* Perform LU decomposition for matrix1 */
+ if (LUDecompose(&matrix1[0][0], 5, 5) < 0) {
+ fprintf(stderr, "Error during LU decomposition for matrix1.\n");
+ return -1;
+ }
+
+ /* Calculate determinant for matrix1 */
+ det1 = LUDeterminant(&matrix1[0][0], 5, 5);
+
+ matrix2[0][0] = ABD[1][1];
+ matrix2[0][1] = ABD[1][2];
+ matrix2[0][2] = ABD[1][3];
+ matrix2[0][3] = ABD[1][4];
+ matrix2[0][4] = ABD[1][5];
+ matrix2[1][0] = ABD[2][1];
+ matrix2[1][1] = ABD[2][2];
+ matrix2[1][2] = ABD[2][3];
+ matrix2[1][3] = ABD[2][4];
+ matrix2[1][4] = ABD[2][5];
+ matrix2[2][0] = ABD[3][1];
+ matrix2[2][1] = ABD[3][2];
+ matrix2[2][2] = ABD[3][3];
+ matrix2[2][3] = ABD[3][4];
+ matrix2[2][4] = ABD[3][5];
+ matrix2[3][0] = ABD[4][1];
+ matrix2[3][1] = ABD[4][2];
+ matrix2[3][2] = ABD[4][3];
+ matrix2[3][3] = ABD[4][4];
+ matrix2[3][4] = ABD[4][5];
+ matrix2[4][0] = ABD[5][1];
+ matrix2[4][1] = ABD[5][2];
+ matrix2[4][2] = ABD[5][3];
+ matrix2[4][3] = ABD[5][4];
+ matrix2[4][4] = ABD[5][5];
+
+ // Perform LU decomposition for matrix2
+ if (LUDecompose(&matrix2[0][0], 5, 5) < 0) {
+ fprintf(stderr, "Error during LU decomposition for matrix2.\n");
+ return -1;
+ }
+
+ // Calculate determinant for matrix2
+ det2 = LUDeterminant(&matrix2[0][0], 5, 5);
+
+ // Check for division by zero
+ if (det2 == 0) {
+ fprintf(stderr, "Error: Division by zero in E_x calculation, determinant of the submatrix is zero.\n");
+ return -1;
+ }
+
+ // Compute E_x
+ nxy = -det1 / det2;
+ return nxy;
+}
+double nu_yx(double ABD[6][6], double h) {
+ double matrix1[5][5], matrix2[5][5];
+ double det1, det2, nyx;
+
+ matrix1[0][0] = ABD[0][1];
+ matrix1[0][1] = ABD[0][2];
+ matrix1[0][2] = ABD[0][3];
+ matrix1[0][3] = ABD[0][4];
+ matrix1[0][4] = ABD[0][5];
+ matrix1[1][0] = ABD[2][0];
+ matrix1[1][1] = ABD[2][2];
+ matrix1[1][2] = ABD[2][3];
+ matrix1[1][3] = ABD[2][4];
+ matrix1[1][4] = ABD[2][5];
+ matrix1[2][0] = ABD[3][0];
+ matrix1[2][1] = ABD[3][2];
+ matrix1[2][2] = ABD[3][3];
+ matrix1[2][3] = ABD[3][4];
+ matrix1[2][4] = ABD[3][5];
+ matrix1[3][0] = ABD[4][0];
+ matrix1[3][1] = ABD[4][2];
+ matrix1[3][2] = ABD[4][3];
+ matrix1[3][3] = ABD[4][4];
+ matrix1[3][4] = ABD[4][5];
+ matrix1[4][0] = ABD[5][0];
+ matrix1[4][1] = ABD[5][2];
+ matrix1[4][2] = ABD[5][3];
+ matrix1[4][3] = ABD[5][4];
+ matrix1[4][4] = ABD[5][5];
+
+ // Perform LU decomposition for matrix1
+ if (LUDecompose(&matrix1[0][0], 5, 5) < 0) {
+ fprintf(stderr, "Error during LU decomposition for matrix1.\n");
+ return -1;
+ }
+
+ // Calculate determinant for matrix1
+ det1 = LUDeterminant(&matrix1[0][0], 5, 5);
+
+ matrix2[0][0] = ABD[0][0];
+ matrix2[0][1] = ABD[0][2];
+ matrix2[0][2] = ABD[0][3];
+ matrix2[0][3] = ABD[0][4];
+ matrix2[0][4] = ABD[0][5];
+ matrix2[1][0] = ABD[2][0];
+ matrix2[1][1] = ABD[2][2];
+ matrix2[1][2] = ABD[2][3];
+ matrix2[1][3] = ABD[2][4];
+ matrix2[1][4] = ABD[2][5];
+ matrix2[2][0] = ABD[3][0];
+ matrix2[2][1] = ABD[3][2];
+ matrix2[2][2] = ABD[3][3];
+ matrix2[2][3] = ABD[3][4];
+ matrix2[2][4] = ABD[3][5];
+ matrix2[3][0] = ABD[4][0];
+ matrix2[3][1] = ABD[4][2];
+ matrix2[3][2] = ABD[4][3];
+ matrix2[3][3] = ABD[4][4];
+ matrix2[3][4] = ABD[4][5];
+ matrix2[4][0] = ABD[5][0];
+ matrix2[4][1] = ABD[5][2];
+ matrix2[4][2] = ABD[5][3];
+ matrix2[4][3] = ABD[5][4];
+ matrix2[4][4] = ABD[5][5];
+
+ // Perform LU decomposition for matrix2
+ if (LUDecompose(&matrix2[0][0], 5, 5) < 0) {
+ fprintf(stderr, "Error during LU decomposition for matrix2.\n");
+ return -1;
+ }
+
+ // Calculate determinant for matrix2
+ det2 = LUDeterminant(&matrix2[0][0], 5, 5);
+
+ // Check for division by zero
+ if (det2 == 0) {
+ fprintf(stderr, "Error: Division by zero in E_x calculation, determinant of the submatrix is zero.\n");
+ return -1;
+ }
+
+ // Compute E_x
+ nyx = -det1 / det2;
+ return nyx;
+}
+
+/*void*/
+/*calc_stressCLT(*/
+/* MaterialProperties *mat_list, Laminate *lam, */
+/* double F[6], int *fail_list, double dT, double dM, */
+/* StressStrainResult *result*/
+/* ) */
+/*{*/
+/* int num = lam->num_layers;*/
+/* double *Z = malloc((num + 1) * sizeof(double));*/
+/* assemble_Z(Z, &lam);*/
+/**/
+/* double h; for(i = 0; i < num; i++) h += lam.thk[i];*/
+/**/
+/* double ABD[N][N];*/
+/* assemble_ABD(ABD, materials, &lam, Z, NULL);*/
+/**/
+/* printf("%d\n", NELEMS(Z));*/
+/*}*/
diff --git a/laminate.dat b/laminate.dat
@@ -0,0 +1,7 @@
+1 0.127 0
+1 0.127 90
+0 0.127 0
+2 0.127 0
+0 0.127 0
+1 0.127 90
+1 0.127 0
diff --git a/main b/main
Binary files differ.
diff --git a/main.c b/main.c
@@ -0,0 +1,157 @@
+/*
+ * main.c
+ *
+ * Author: Julio C. Ramirez-Ceballos
+ * Licence: MIT
+ * Description:
+ * This code is used to calculate the equivalent properties of a laminate
+ * with some configuration.
+ * INPUT:
+ * - Material properties
+ * - Layer
+ */
+
+#include <stdio.h>
+#include "clt.h"
+
+int main()
+{
+ int i;
+
+ /*
+ * ADD A FUNCTION THAT DETERMINES THIS VALUES FROM THE MICRO-SCALE
+ * (FIBER + MATRIX & VOIDS)
+ */
+
+ double F[] = {
+ 100000,
+ 0e2,
+ 0e4,
+ 10e0,
+ 0e0,
+ 0e0,
+ };
+
+ int num_materials;
+ MaterialProperties *materials =
+ read_material_properties("materials.dat", &num_materials, true);
+
+ Laminate lam;
+ read_laminate_data("laminate.dat", &lam, true);
+
+ double *Z = malloc((lam.num_layers + 1) * sizeof(double));
+
+ assemble_Z(Z, &lam);
+ printf("Layer Interface Positions (Z):\n");
+ for (i = 0; i <= lam.num_layers; ++i) {
+ printf("Z[%d] = %.4e mm\n", i, Z[i]);
+ }
+
+ printf("\n");
+
+ double ABD[N][N];
+ assemble_ABD(ABD, materials, &lam, Z, NULL);
+
+ if (LUDecompose(&ABD[0][0], 6, 6) < 0) {
+ fprintf(stderr, "Error during LU decomposition for matrix1.\n");
+ return -1;
+ }
+
+ LUSolve(&ABD[0][0], &F[0], 6, 6);
+
+ printf("\n");
+ printf("Strain vector:\n");
+ for (i = 0; i < 6; i++) {
+ printf("%e\n", F[i]);
+ }
+
+ double *strain = malloc(3 * sizeof(double));
+ double *curvature = malloc(3 * sizeof(double));
+
+ for (i = 0; i < 3; i++) strain[i] = F[i];
+ for (i = 0; i < 3; i++) curvature[i] = F[i+3];
+
+ /* Allocate and initialize Laminate System strain vectors */
+ double LS_strain_sup[3][lam.num_layers];
+ double LS_strain_inf[3][lam.num_layers];
+
+ for (int i = 0; i < lam.num_layers; i++) {
+ for(int j = 0; j < 3; j++) {
+ LS_strain_inf[j][i] = strain[j] + curvature[j] * Z[i];
+ LS_strain_sup[j][i] = strain[j] + curvature[j] * Z[i + 1];
+ }
+ }
+ printf("Laminate System Strains:\n");
+ printf("Layer | Inferior (Z-) | Superior (Z+)\n");
+ printf("------|----------------------------------|----------------------------\n");
+ for (int i = 0; i < lam.num_layers; i++) {
+ printf("%4d | ", i + 1);
+ for (int j = 0; j < 3; j++) {
+ printf("%.3e ", LS_strain_inf[j][i]);
+ }
+ printf(" | ");
+ for (int j = 0; j < 3; j++) {
+ printf("%.3e ", LS_strain_sup[j][i]);
+ }
+ printf("\n");
+ }
+ printf("\n");
+
+ double MS_strain_sup[3][lam.num_layers];
+ double MS_strain_inf[3][lam.num_layers];
+ double MS_stress_sup[3][lam.num_layers];
+ double MS_stress_inf[3][lam.num_layers];
+ double T[3][3], Q[3][3];
+
+
+
+ /*CORRIGER*/
+ /*for (i = 0; i < lam.num_layers; i++) {*/
+ /* MaterialProperties *mat_prop = &materials[lam.mat_id[i]];*/
+ /* assemble_Q(Q, mat_prop, fail_list ? fail_list[i] : 0);*/
+ /* assemble_T(T, lam.ang[i]);*/
+ /* double temp[3];*/
+ /**/
+ /* dot_product(temp, T, LS_strain_sup[0][i], 3, 3, 1);*/
+ /* for(int j = 0; j < 3; j++) MS_strain_sup[j][i] = temp[j];*/
+ /**/
+ /* dot_product(temp, T, LS_strain_inf[0][i], 3, 3, 1);*/
+ /* for(int j = 0; j < 3; j++) MS_strain_sup[j][i] = temp[j];*/
+ /**/
+ /* dot_product(temp, Q, MS_strain_sup[0][i], 3, 3, 1);*/
+ /* for(int j = 0; j < 3; j++) MS_stress_sup[j][i] = temp[j];*/
+ /**/
+ /* dot_product(temp, Q, MS_strain_inf[0][i], 3, 3, 1);*/
+ /* for(int j = 0; j < 3; j++) MS_stress_inf[j][i] = temp[j];*/
+ /*}*/
+ /*printf("%.4e\n", MS_stress_inf[0][0]);*/
+
+ double h;
+ for(i = 0; i < lam.num_layers; i++) {
+ h += lam.thk[i];
+ }
+ printf("h = %0.4e mm\n", h);
+
+ double Ex = E_x(ABD, h);
+ printf("Equivalent elastic modulus Ex:\t%.4e\n", Ex);
+
+ double Ey = E_y(ABD, h);
+ printf("Equivalent elastic modulus Ey:\t%.4e\n", Ey);
+
+ double Gxy = G_xy(ABD, h);
+ printf("Equivalent shear modulus Gxy:\t%.4e\n", Gxy);
+
+ double nxy = nu_xy(ABD, h);
+ printf("Equivalent poisson ratio nu_xy:\t%.2f\n", nxy);
+
+ double nyx = nu_yx(ABD, h);
+ printf("Equivalent poisson ratio nu_yx:\t%.2f\n", nyx);
+
+ free(lam.thk);
+ free(lam.ang);
+ free(lam.mat_id);
+ free(Z);
+
+ return 0;
+}
+
diff --git a/materials.dat b/materials.dat
@@ -0,0 +1,14 @@
+# 'I' = ISOTROPIC
+# I E nu
+
+# 'A' = ANISOTROPIC
+# A E1 E2 nu12 G12 a1 a2 b1 b2 Xt Xc Yt Yc S12 S23
+
+# Rubber
+I 69e6 0.354
+
+# Carbon
+A 69e6 6e6 0.354 3e6 2.1e-3 2.1e-3 0.01 0.35 47e3 14e3 24e3 18e3 75e3 41e3 2e4
+
+# Steel
+I 210000e6 0.354
diff --git a/rw.c b/rw.c
@@ -0,0 +1,218 @@
+/*
+ * rw.c
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdbool.h>
+
+#define MAX_mat 100 // Assuming a maximum number of mat
+#define MAX_LINE_LENGTH 256 // Assuming lines won't exceed this length
+
+/* Structure for material properties */
+typedef struct {
+ double E1, E2, nu12, G12,
+ a1, a2, b1, b2,
+ Xt, Xc, Yt, Yc, S12, S32;
+ bool isotropic;
+} MaterialProperties;
+
+/* Structure for laminate properties */
+typedef struct {
+ int *mat_id;
+ double *thk;
+ double *ang;
+ int num_layers;
+} Laminate;
+
+void read_laminate_data(const char *filename, Laminate *lam, bool print_details) {
+ FILE *file;
+ char line[100]; // Assuming each line won't exceed 100 characters
+ int layer_count = 0;
+
+ // Attempt to open the file
+ file = fopen(filename, "r");
+ if (file == NULL) {
+ fprintf(stderr, "Error opening file %s\n", filename);
+ exit(1);
+ }
+
+ // Count the number of lines to determine num_layers
+ while (fgets(line, sizeof(line), file) != NULL) {
+ if (line[0] != '\n' && line[0] != '#') { // Ignore empty lines or comments
+ layer_count++;
+ }
+ }
+
+ // Allocate memory for the laminate structure
+ lam->num_layers = layer_count;
+ lam->mat_id = malloc(layer_count * sizeof(int));
+ lam->thk = malloc(layer_count * sizeof(double));
+ lam->ang = malloc(layer_count * sizeof(double));
+
+ if (lam->mat_id == NULL || lam->thk == NULL || lam->ang == NULL) {
+ fprintf(stderr, "Memory allocation failed\n");
+ fclose(file);
+ exit(1);
+ }
+
+ // Reset file pointer to beginning of file
+ rewind(file);
+
+ // Read the actual data
+ int i = 0;
+ while (fgets(line, sizeof(line), file) != NULL) {
+ if (line[0] != '\n' && line[0] != '#') { // Skip empty lines or comments
+ // Temporary variables to hold the parsed data from each line
+ int mat_id;
+ double thickness, angle;
+
+ // Use sscanf to parse the line
+ if (sscanf(line, "%d %lf %lf", &mat_id, &thickness, &angle) == 3) {
+ lam->mat_id[i] = mat_id;
+ lam->thk[i] = thickness;
+ lam->ang[i] = angle;
+ i++;
+ } else {
+ fprintf(stderr, "Error reading line: %s\n", line);
+ // Optionally, you could choose to continue or break here
+ }
+ }
+ }
+
+ // Print the laminate details if requested
+ if (print_details) {
+ printf("\n");
+ printf("+------------------------+\n");
+ printf("| Laminate Configuration |\n");
+ printf("+------------------------+\n\n");
+ printf("Number of Layers: %d\n", lam->num_layers);
+ printf("Layer | Thickness (mm) | Angle (°) | Material ID\n");
+ printf("------|----------------|-----------|------------\n");
+ for (int i = 0; i < lam->num_layers; i++) {
+ printf(" %2d | %13.3f | %9.2f | %10d\n",
+ i + 1,
+ lam->thk[i],
+ lam->ang[i],
+ lam->mat_id[i]);
+ }
+ printf("\n");
+ }
+ fclose(file);
+}
+
+MaterialProperties*
+read_material_properties(const char *filename,
+ int *n_mat,
+ bool print_details)
+{
+ FILE *file;
+ char line[MAX_LINE_LENGTH];
+ int material_count = 0;
+ MaterialProperties *mat = NULL;
+
+ printf("\n");
+ printf("+---------------------+\n");
+ printf("| Material Properties |\n");
+ printf("+---------------------+\n");
+ printf("\n");
+
+ file = fopen(filename, "r");
+ if (file == NULL) {
+ fprintf(stderr, "Error opening file %s\n", filename);
+ return NULL;
+ }
+
+ // First pass: count the number of mat
+ while (fgets(line, sizeof(line), file)) {
+ if (line[0] == '#' || line[0] == '\n' || line[0] == '\r') continue;
+ material_count++;
+ }
+
+ // Reset file pointer to beginning
+ rewind(file);
+
+ // Allocate memory for mat
+ *n_mat = material_count;
+ mat = (MaterialProperties *)malloc(
+ material_count * sizeof(MaterialProperties));
+ if (mat == NULL) {
+ fprintf(stderr, "Memory allocation failed\n");
+ fclose(file);
+ return NULL;
+ }
+
+ // Second pass: read the material properties
+ int mat_indx = 0;
+ while (fgets(line, sizeof(line), file)) {
+ if (line[0] == '#' || line[0] == '\n' || line[0] == '\r') continue;
+
+ char type;
+ // Check if the material type is specified at the beginning of the line
+ if (sscanf(line, "%c", &type) != 1) {
+ fprintf(stderr, "Error reading material type from line: %s", line);
+ continue;
+ }
+
+ if (type == 'I') { // 'I' for Isotropic
+ // For isotropic mat, we only need E and nu (Poisson's ratio)
+ if (sscanf(line + 1, "%lf %lf",
+ &mat[mat_indx].E1,
+ &mat[mat_indx].nu12) != 2) {
+ fprintf(stderr, "Error reading isotropic material properties from line: %s", line);
+ continue;
+ }
+ // Set isotropic flag
+ mat[mat_indx].isotropic = true;
+ // Fill in other properties with default or derived values
+ mat[mat_indx].E2 = mat[mat_indx].E1; // For isotropic mat, E1 = E2
+ mat[mat_indx].G12 = mat[mat_indx].E1 / (2.0 * (1 + mat[mat_indx].nu12));
+ if (print_details) {
+ printf("ISOTROPIC:\n");
+ printf("Material ID %d Properties:\n", mat_indx);
+ printf(" E (Pa) : %.2e\n", mat[mat_indx].E1);
+ printf(" nu : %.3f\n\n", mat[mat_indx].nu12);
+ }
+ } else if (type == 'A') { // 'A' for Anisotropic
+ // Read all properties for anisotropic mat
+ if (sscanf(line + 1, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
+ &mat[mat_indx].E1, &mat[mat_indx].E2,
+ &mat[mat_indx].nu12, &mat[mat_indx].G12,
+ &mat[mat_indx].a1, &mat[mat_indx].a2,
+ &mat[mat_indx].b1, &mat[mat_indx].b2,
+ &mat[mat_indx].Xt, &mat[mat_indx].Xc,
+ &mat[mat_indx].Yt, &mat[mat_indx].Yc,
+ &mat[mat_indx].S12, &mat[mat_indx].S32) != 14) {
+ fprintf(stderr, "Error reading anisotropic material properties from line: %s", line);
+ continue;
+ }
+ // Set isotropic flag to false
+ mat[mat_indx].isotropic = false;
+
+ if (print_details) {
+ printf("ORTHOTROPIC:\n");
+ printf("Material ID %d Properties:\n", mat_indx);
+ printf(" E1 (Pa) : %.2e\n", mat[mat_indx].E1);
+ printf(" E2 (Pa) : %.2e\n", mat[mat_indx].E2);
+ printf(" nu12 : %.3f\n", mat[mat_indx].nu12);
+ printf(" G12 (Pa) : %.2e\n", mat[mat_indx].G12);
+ printf(" a1 (1/K) : %.2e\n", mat[mat_indx].a1);
+ printf(" a2 (1/K) : %.2e\n", mat[mat_indx].a2);
+ printf(" b1 : %.2e\n", mat[mat_indx].b1);
+ printf(" b2 : %.2e\n", mat[mat_indx].b2);
+ printf(" Xt (Pa) : %.2e\n", mat[mat_indx].Xt);
+ printf(" Xc (Pa) : %.2e\n", mat[mat_indx].Xc);
+ printf(" Yt (Pa) : %.2e\n", mat[mat_indx].Yt);
+ printf(" Yc (Pa) : %.2e\n", mat[mat_indx].Yc);
+ printf(" S12 (Pa) : %.2e\n", mat[mat_indx].S12);
+ printf(" S23 (Pa) : %.2e\n\n", mat[mat_indx].S32);
+ }
+ } else {
+ fprintf(stderr, "Unknown material type in line: %s", line);
+ continue;
+ }
+ mat_indx++;
+ }
+ fclose(file);
+ return mat;
+}
diff --git a/utils.h b/utils.h
@@ -0,0 +1,260 @@
+/*
+ * utils.h
+ */
+
+#include <math.h>
+
+#define PI 3.14159265358979323846
+#define DEGTORAD 0.017453292519943295 /* pi/180 */
+#define CONST13 0.333333333333333
+#define MAT_ELEM(mat, row, col) mat[row][col]
+#define SMALLEST_PIVOT 1.0e-5 /* Smallest allowed non-singular pivot. */
+#define SIZE 6
+
+static int *pivot = NULL; /* Pivot vector. */
+
+/* Simple factorial function for Taylor series */
+double factorial(int n) {
+ if (n == 0 || n == 1) return 1;
+ return n * factorial(n - 1);
+}
+
+/* Sine function using Taylor series */
+double Sin(double x) {
+ double term, sum = 0;
+ int n;
+ x = fmod(x, 2 * PI); /* Reduce x to be within[0, 2π) */
+
+ for (n = 0; n < 10; ++n) { /* Adjust the number of terms for precision */
+ term = pow(-1, n) * pow(x, 2*n+1) / factorial(2*n+1);
+ sum += term;
+ }
+ return sum;
+}
+
+/* Cosine function using Taylor series */
+double Cos(double x) {
+ double term, sum = 0;
+ int n;
+ x = fmod(x, 2 * PI); /* Reduce x to be within[0, 2π) */
+
+ /* Adjust the number of terms for precision */
+ for (n = 0; n < 10; ++n) {
+ term = pow(-1, n) * pow(x, 2*n) / factorial(2*n);
+ sum += term;
+ }
+ return sum;
+}
+
+void matmul(double A[3][3], double B[3][3], double result[3][3])
+{
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ result[i][j] = 0.0;
+ for (int k = 0; k < 3; k++) {
+ result[i][j] += A[i][k] * B[k][j];
+ }
+ }
+ }
+}
+
+/*CORRIGER*/
+/*void dot_product(double result[3][1], double A[3][3], double B[3][1], int m, int n, int p) {*/
+/* for(int i = 0; i < m; ++i) {*/
+/* for(int j = 0; j < p; ++j) {*/
+/* result[i*p + j] = 0.0;*/
+/* for(int k = 0; k < n; ++k) {*/
+/* result[i*p + j] += A[i*n + k] * B[k*p + j];*/
+/* }*/
+/* }*/
+/* }*/
+/*}*/
+
+void transpose(double A[3][3], double Ai[3][3])
+{
+ for (int j = 0; j < 3; j++) {
+ for (int k = 0; k < 3; k++) {
+ Ai[j][k] = A[k][j];
+ }
+ }
+}
+
+ void
+dispmat(double *A, int rows, int cols)
+{
+ for (int i = 0; i < rows; i++) {
+ for (int j = 0; j < cols; j++) {
+ /*printf("%.2e\t", A[i][j]);*/
+ printf("%.2e\t", *(A + i * cols + j));
+ }
+ printf("\n");
+ }
+}
+
+/*
+ * This file contains routines for solving an nxn system of linear equations
+ * using LU decomposition. The implementation is based on material in
+ * "Numerical Methods - A Software Approach" by Johnston.
+ * Function LUDecompose performs the decomposition step, LUSolve solves a
+ * system given a particular right hand side vector b, and LUDetermiant will
+ * return the determinant of a matrix which has been decomposed to LU form.
+ */
+
+int LUDecompose(double *amat, int n, int numcols)
+{
+
+ /*
+ * This function performs the LU decomposition step. The coefficient
+ * matrix 'amat' will be overwritten with the LU matrices in the
+ * process.
+ *
+ * Parameters :
+ * amat - This is the coefficient matrix, where the coefficients
+ * are doubles. This function performs the row-major mapping
+ * itself, declaring 'amat' to be a pointer rather than an array.
+ * The calling routine can define the actual parameter to be
+ * either.
+ *
+ * n - This parameter indicates the size of the current system
+ * (nxn).
+ *
+ * numcols - This parameter indicates the actual defined column
+ * dimension of the coefficient matrix. A bit reminiscent of
+ * FORTRAN.
+ */
+
+ int i, j, k, n_minus_1 = n - 1;
+ double dtmp1, *dptr1, *dptr2;
+
+ /* To keep track of determinant sign due to row swaps*/
+ int sign = 1;
+
+ if (pivot != NULL) free(pivot);
+ if ((pivot = (int *) malloc(n * sizeof(int))) == NULL) {
+ fprintf(stderr, "Error in LUDecompose - malloc \n");
+ return -2;
+ }
+
+ /* Initialize pivot vector*/
+ for (i = 0; i < n_minus_1; i++) {
+ *(pivot + i) = i;
+ }
+ *(pivot + n_minus_1) = n_minus_1;
+
+ for (i = 0; i < n_minus_1; i++) {
+ k = i;
+ dptr1 = amat + i * numcols + i;
+ for (j = i + 1; j < n; j++) {
+ dptr2 = amat + j * numcols + i;
+ if (fabs(*dptr2) > fabs(*dptr1) ) {
+ dptr1 = dptr2;
+ k = j;
+ }
+ }
+
+ if (fabs(*dptr1) < SMALLEST_PIVOT) {
+ fprintf(stderr, "Error in LUDecompose - matrix is singular\n");
+ return -1;
+ }
+
+ /* Swap rows and update sign */
+ if (k != i) {
+ int temp = *(pivot + i);
+ *(pivot + i) = *(pivot + k);
+ *(pivot + k) = temp;
+
+ /* Flip the sign of the determinant for each swap */
+ sign *= -1;
+
+ for (j = 0; j < n; j++) {
+ dtmp1 = *(amat + i * numcols + j);
+ *(amat + i * numcols + j) = *(amat + k * numcols + j);
+ *(amat + k * numcols + j) = dtmp1;
+ }
+ }
+
+ for (j = i+1; j < n; j++) {
+ dtmp1 = *(amat + j * numcols + i) / *dptr1;
+ *(amat + j * numcols + i) = dtmp1;
+ for (k = i+1; k < n; k++) {
+ *(amat + j * numcols + k) -= dtmp1 * *(amat + i * numcols + k);
+ }
+ }
+ }
+ *(pivot + n - 1) = sign;
+ return 1;
+}
+
+void LUSolve(amat, b, n, numcols)
+ double *amat, *b;
+ int n, numcols;
+ /*
+ * This function solves a system of equations given a particular right
+ * hand side vector 'b'. The coefficient matrix 'amat' must be
+ * decomposed by the function LUDecompose prior to calling LUSolve.
+ *
+ * Parameters:
+ * amat - This is the LU decomposition of the coefficient matrix.
+ * The original coefficient matrix must first be passed to
+ * the function LUDecompose before calling this routine.
+ *
+ * b - This parameter is the right hand side vector b. This
+ * routine will compute the solution vector and return this
+ * result in the parameter 'b'. NOTE - b WILL BE OVERWRITTEN
+ *
+ * n - This parameter indicates the size of the current system
+ * (nxn).
+ *
+ * numcols - This parameter indicates the actual defined column
+ * dimension of the coefficient matrix. A bit
+ * reminiscent of FORTRAN.
+ */
+{
+ double dtmp1, *dptr1;
+ int i, j, n_minus_1 = n - 1;
+
+ /* Perform the row interchanges recorded in 'pivot' to vector b. */
+ for (i = 0; i < n; i++) {
+ /* Interchange element i and j. */
+ j = *(pivot + i);
+ if (j != i) {
+ dtmp1 = *(b + i);
+ *(b + i) = *(b + j);
+ *(b + j) = dtmp1;
+ }
+ }
+
+ /* Solve the system Ld = Pb for d. Vector d will overwritte b. */
+ for (i = 0; i < n; i++) {
+ dtmp1 = *(b + i);
+ for (j = 0, dptr1 = amat + i*numcols; j < i; j++, dptr1++)
+ dtmp1 -= *dptr1 * *(b + j);
+ *(b + i) = dtmp1;
+ }
+
+ /* Solve the system Ux = d for x. Vector x will overwritte b (and d). */
+ for (i = n_minus_1; i >= 0; i--) {
+ dtmp1 = *(b + i);
+ dptr1 = amat + i * numcols + n_minus_1;
+
+ for (j = n_minus_1; j > i; j--) {
+ dtmp1 -= *dptr1 * *(b + j);
+ dptr1--;
+ }
+ *(b + i) = dtmp1 / *dptr1;
+ }
+}
+
+double LUDeterminant(double *amat, int n, int numcols)
+ /*
+ * This function will compute the determinant of the matrix 'amat'.
+ * 'amat' must have already been decomposed to LU form (use
+ * LUDecompose).
+ */
+{
+ double det = 1.0;
+ for (int i = 0; i < n; i++) {
+ det *= *(amat + i * (numcols + 1));
+ }
+ return det * *(pivot + n - 1);
+}
