20 using namespace Hermes::Algebra::DenseMatrixOperations;
26 double* Shapeset::calculate_constrained_edge_combination(
int order,
int part,
int ori, ElementMode2D mode)
33 assert((order - ebias) >= 0);
39 for (n = 2; n <= part; n <<= 1)
43 double hi = -((double) part * n2 - 1.0);
44 double lo = -((double) (part + 1) * n2 - 1.0);
48 for (i = 0; i <= order; i++)
49 idx[i] = get_edge_index(0, ori, i, mode);
53 double f_lo = 0.0, f_hi = 0.0;
56 f_lo = get_value(0, idx[order], lo, -1.0, 0, mode);
57 f_hi = get_value(0, idx[order], hi, -1.0, 0, mode);
67 n = order + 1 - ebias;
68 int space_type = this->get_space_type();
69 int component = (space_type == HERMES_HDIV_SPACE)? 1 : 0;
70 double** a = new_matrix<double>(n, n);
71 double* b =
new double[n];
72 for (i = 0; i < n; i++)
75 int o = (ebias == 0) ? order + 1 : order;
76 double p = cos((i + 1) * M_PI / o);
77 double r = (p + 1.0) * 0.5;
81 for (j = 0; j < n; j++)
82 a[i][j] = get_value(0, idx[j + ebias], p, -1.0, component, mode);
85 b[i] = c * get_value(0, idx[order], lo*s + hi*r, -1.0, component, mode) - f_lo*s - f_hi*r;
90 int* iperm =
new int[n];
91 ludcmp(a, n, iperm, &d);
92 lubksb(a, n, iperm, b);
101 double* Shapeset::get_constrained_edge_combination(
int order,
int part,
int ori,
int& nitems, ElementMode2D mode)
103 int index = 2*((max_order + 1 - ebias)*part + (order - ebias)) + ori;
106 if(comb_table == NULL)
109 while (table_size <= index) table_size *= 2;
110 comb_table = (
double**) malloc(table_size *
sizeof(
double*));
111 memset(comb_table, 0, table_size *
sizeof(
double*));
113 else if(index >= table_size)
116 int old_size = table_size;
117 while (index >= table_size) table_size *= 2;
120 comb_table = (
double**) realloc(comb_table, table_size *
sizeof(
double*));
121 memset(comb_table + old_size, 0, (table_size - old_size) *
sizeof(
double*));
125 if(comb_table[index] == NULL)
128 comb_table[index] = calculate_constrained_edge_combination(order, part, ori, mode);
131 nitems = order + 1 - ebias;
132 return comb_table[index];
135 void Shapeset::free_constrained_edge_combinations()
137 if(comb_table != NULL)
139 for (
int i = 0; i < table_size; i++)
140 if(comb_table[i] != NULL)
141 delete [] comb_table[i];
148 double Shapeset::get_constrained_value(
int n,
int index,
double x,
double y,
int component, ElementMode2D mode)
152 int part = (unsigned) index >> 7;
153 int order = (index >> 3) & 15;
154 int edge = (index >> 1) & 3;
158 double sum, *comb = get_constrained_edge_combination(order, part, ori, nc, mode);
161 shape_fn_t* table = shape_table[n][mode][component];
162 for (i = 0; i < nc; i++)
163 sum += comb[i] * table[get_edge_index(edge, ori, i + ebias, mode)](x, y);
168 Shapeset::~Shapeset() { free_constrained_edge_combinations(); }
170 int Shapeset::get_max_order()
const {
return max_order; }
172 int Shapeset::get_num_components()
const {
return num_components; }
174 int Shapeset::get_vertex_index(
int vertex, ElementMode2D mode)
const
176 if(mode == HERMES_MODE_TRIANGLE)
180 return vertex_indices[mode][vertex];
183 int Shapeset::get_edge_index(
int edge,
int ori,
int order, ElementMode2D mode)
const
185 if(mode == HERMES_MODE_TRIANGLE)
189 assert(order >= 0 && order <= max_order);
190 assert(ori == 0 || ori == 1);
191 return edge_indices[mode][edge][2*order + ori];
194 int* Shapeset::get_bubble_indices(
int order, ElementMode2D mode)
const
197 assert(H2D_GET_V_ORDER(order) >= 0 && H2D_GET_V_ORDER(order) <= max_order);
199 if(mode == HERMES_MODE_QUAD)
200 index = H2D_MAKE_QUAD_ORDER(H2D_GET_V_ORDER(order),
H2D_GET_H_ORDER(order));
201 return bubble_indices[mode][index];
204 int Shapeset::get_num_bubbles(
int order, ElementMode2D mode)
const
207 assert(H2D_GET_V_ORDER(order) >= 0 && H2D_GET_V_ORDER(order) <= max_order);
208 return bubble_count[mode][order];
211 int Shapeset::get_constrained_edge_index(
int edge,
int order,
int ori,
int part, ElementMode2D mode)
const
213 if(mode == HERMES_MODE_TRIANGLE)
217 assert(order >= 0 && order <= max_order);
219 assert(order <= H2D_ORDER_MASK);
220 return -1 - ((part << 7) + (order << 3) + (edge << 1) + ori);
223 int Shapeset::get_order(
int index, ElementMode2D mode)
const
226 return index_to_order[mode][index];
228 else return ((-1 - index) >> 3) & 15;
231 double Shapeset::get_value(
int n,
int index,
double x,
double y,
int component, ElementMode2D mode)
235 Shapeset::shape_fn_t** shape_expansion = shape_table[n][mode];
236 if(shape_expansion == NULL)
239 static int warned_mode = -1, warned_index = -1, warned_n = 1;
240 this->warn_if(warned_mode != mode || warned_index != index || warned_n != n,
"Requested undefined expansion %d (mode: %d) of a shape %d, returning 0", n, mode, index);
242 warned_index = index;
247 return shape_expansion[component][index](x, y);
250 return get_constrained_value(n, index, x, y, component, mode);
253 double Shapeset::get_fn_value (
int index,
double x,
double y,
int component, ElementMode2D mode) {
return get_value(0, index, x, y, component, mode); }
254 double Shapeset::get_dx_value (
int index,
double x,
double y,
int component, ElementMode2D mode) {
return get_value(1, index, x, y, component, mode); }
255 double Shapeset::get_dy_value (
int index,
double x,
double y,
int component, ElementMode2D mode) {
return get_value(2, index, x, y, component, mode); }
256 double Shapeset::get_dxx_value(
int index,
double x,
double y,
int component, ElementMode2D mode) {
return get_value(3, index, x, y, component, mode); }
257 double Shapeset::get_dyy_value(
int index,
double x,
double y,
int component, ElementMode2D mode) {
return get_value(4, index, x, y, component, mode); }
258 double Shapeset::get_dxy_value(
int index,
double x,
double y,
int component, ElementMode2D mode) {
return get_value(5, index, x, y, component, mode); }
261 double2* Shapeset::get_ref_vertex(
int vertex, ElementMode2D mode)
263 return &ref_vert[mode][vertex];