linear_system_full.c 5.04 KB
 Jonathan Lambrechts committed Feb 02, 2017 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 ``````#include #include "mesh.h" #include "linear_system.h" #include #include //from wikipedia //INPUT: A - array of pointers to rows of a square matrix having dimension N // Tol - small tolerance number to detect failure when the matrix is near degenerate //OUTPUT: Matrix A is changed, it contains both matrices L-E and U as A=(L-E)+U such that P*A=L*U. // P - array of N+1 integers containing pivoting of A and P[N] is for determinant computation // (P[i] is an index of the row in A placed at the i-th row static int LUPDecompose(double **A,int N,double Tol,int *P){ int i,j,k,imax; double maxA,*ptr,absA; for(i=0;imaxA){ maxA=absA; imax=k; } if(maxA=0;i--){ for(int k=i+1;kmesh; system->n = mesh->n_nodes*system->n_fields; system->isfixed = malloc(sizeof(int)*system->n); for (int i = 0; i < system->n; ++i) { system->isfixed[i] = 0; } for (int ibnd = 0; ibnd < n_boundaries; ++ibnd) { const StrongBoundary *bnd = bnds + ibnd; int iphys; for (iphys = 0; iphys < mesh->n_phys; ++iphys) { if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0) break; } if (iphys == mesh->n_phys) printf("Boundary tag \"%s\" not found.", bnd->tag); for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){ system->isfixed[mesh->phys_nodes[iphys][i]*system->n_fields+bnd->field] = 1; } } } LinearSystem *linear_system_new(Mesh *mesh, int n_fields, int n_boundaries, const StrongBoundary *boundaries) { LinearSystem *system = malloc(sizeof(LinearSystem)); system->n_fields = n_fields; system->mesh = mesh; linear_system_create_map(system, n_boundaries, boundaries); system->A = malloc(sizeof(double)*system->n*system->n); system->b = malloc(sizeof(double)*system->n); system->x = malloc(sizeof(double)*system->n); return system; } void linear_system_add_to_matrix(LinearSystem *system, int el0, int el1, const double *local_matrix){ int *tri0 = &system->mesh->triangles[el0*3]; int *tri1 = &system->mesh->triangles[el1*3]; int nf = system->n_fields; for (int i = 0; i < 3; ++i) { for (int inf = 0; inf < 3; ++inf) { int ii = tri0[i]*nf + inf; for (int j = 0; j < 3; ++j) { for (int jnf = 0; jnf < nf; ++jnf) { int jj = tri1[j]*nf + jnf; system->A[ii*system->n+jj] += local_matrix[(inf*3+i)*nf*3+jnf*3+j]; } } } } } void linear_system_add_to_rhs(LinearSystem *system, int el0, const double *local_vector) { const int *tri = system->mesh->triangles + el0*3; int n_fields = system->n_fields; for (int i = 0; i < n_fields; ++i) { for (int j = 0; j < 3; ++j) { int m = tri[j]*n_fields+i; system->b[m] += local_vector[i*3+j]; } } } void linear_system_zero_matrix_and_rhs(LinearSystem *system) { for (int i = 0; i < system->n; ++i){ system->b[i] = 0.; for (int j = 0; j < system->n; ++j) { system->A[i*system->n+j] = 0.; } } } void linear_system_solve(LinearSystem *system, double *solution){ int *P = malloc(sizeof(int)*(system->n)); double **rows = malloc(sizeof(double*)*system->n); for (int i = 0; i < system->n; ++i) rows[i] = system->A +i*system->n; for (int i = 0; i < system->n; ++i){ if (system->isfixed[i]) { for (int j = 0; j < system->n; ++j) rows[i][j] = 0.; rows[i][i] = 1.; system->b[i] = 0; } } LUPDecompose(rows, system->n, 1e-8, P); LUPSolve(rows, P, system->b, system->n, system->x); for (int i = 0; i < system->n; ++i) solution[i] = system->x[i]; free(rows); } void linear_system_free(LinearSystem *system) { free(system->b); free(system->x); free(system->A); free(system->isfixed); free(system); } double linear_system_get_rhs_norm(LinearSystem *system) { double norm = 0; for (int i = 0; i < system->n;++i) if (!system->isfixed[i]) norm += system->b[i]*system->b[i]; return sqrt(norm); } void initialize_linear_solver(int argc, char **argv){}``````