69 assert(mp_aff != 0x0) ;
75 if (
etat == ETATZERO) {
84 bool ced = (mgrid.
get_type_r(nzm1) == UNSURR) ;
86 int domain_bc = par_bc.
get_int(0) ;
87 bool bc_ced = ((ced) && (domain_bc == nzm1)) ;
89 int n_conditions = par_bc.
get_int(1) ;
90 assert ((n_conditions==1)||(n_conditions==2)) ;
91 bool derivative = (n_conditions == 2) ;
92 int dl = 0 ;
int l_min = 0 ;
int excised_i =0;
100 bool isexcised = (excised_i==1);
102 int nt = mgrid.
get_nt(0) ;
103 int np = mgrid.
get_np(0) ;
114 int system01_size =0;
116 if (isexcised ==
false){
124 for (
int lz=1; lz<=domain_bc; lz++) {
125 system01_size += n_conditions ;
126 system_size += n_conditions ;
128 assert (
etat != ETATNONDEF) ;
130 bool need_matrices = true ;
133 need_matrices = false ;
135 Matrice system_l0(system01_size, system01_size) ;
136 Matrice system_l1(system01_size, system01_size) ;
137 Matrice system_even(system_size, system_size) ;
138 Matrice system_odd(system_size, system_size) ;
145 int index = 0 ;
int index01 = 0 ;
148 if (isexcised ==
false){
149 system_even.
set(index, index) = 1. ;
150 system_even.
set(index, index + 1) = -1. ;
151 system_odd.
set(index, index) = -(2.*nr - 5.)/alpha ;
152 system_odd.
set(index, index+1) = (2.*nr - 3.)/alpha ;
154 if (domain_bc == 0) {
155 system_l0.
set(index01, index01) = c_val + c_der*4.*(nr-1)*(nr-1)/alpha ;
156 system_l1.
set(index01, index01) = c_val + c_der*(2*nr-3)*(2*nr-3)/alpha ;
158 system_even.
set(index, index-1) = c_val + c_der*4.*(nr-2)*(nr-2)/alpha ;
159 system_even.
set(index, index) = c_val + c_der*4.*(nr-1)*(nr-1)/alpha ;
160 system_odd.
set(index, index-1) = c_val + c_der*(2*nr-5)*(2*nr-5)/alpha ;
161 system_odd.
set(index, index) = c_val + c_der*(2*nr-3)*(2*nr-3)/alpha ;
165 system_l0.
set(index01, index01) = 1. ;
166 system_l1.
set(index01, index01) = 1. ;
167 system_even.
set(index, index-1) = 1. ;
168 system_even.
set(index, index) = 1. ;
169 system_odd.
set(index, index-1) = 1. ;
170 system_odd.
set(index, index) = 1. ;
172 system_l0.
set(index01+1, index01) = 4*(nr-1)*(nr-1)/alpha ;
173 system_l1.
set(index01+1, index01) = (2*nr-3)*(2*nr-3)/alpha ;
174 system_even.
set(index+1, index-1) = 4*(nr-2)*(nr-2)/alpha ;
175 system_even.
set(index+1, index) = 4*(nr-1)*(nr-1)/alpha ;
176 system_odd.
set(index+1, index-1) = (2*nr-5)*(2*nr-5)/alpha ;
177 system_odd.
set(index+1, index) = (2*nr-3)*(2*nr-3)/alpha ;
184 if ((1 == domain_bc)&&(bc_ced))
185 alpha = -0.25/alpha ;
186 system_l0.
set(index01, index01+1) = 1. ;
187 system_l0.
set(index01, index01+2) = -1. ;
188 system_l1.
set(index01, index01+1) = 1. ;
189 system_l1.
set(index01, index01+2) = -1. ;
191 system_even.
set(index, index+1) = 1. ;
192 system_even.
set(index, index+2) = -1. ;
193 system_odd.
set(index, index+1) = 1. ;
194 system_odd.
set(index, index+2) = -1. ;
197 system_l0.
set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
198 system_l0.
set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
199 system_l1.
set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
200 system_l1.
set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
202 system_even.
set(index, index) = -(nr-2)*(nr-2)/alpha ;
203 system_even.
set(index, index+1) = (nr-1)*(nr-1)/alpha ;
204 system_odd.
set(index, index) = -(nr-2)*(nr-2)/alpha ;
205 system_odd.
set(index, index+1) = (nr-1)*(nr-1)/alpha ;
209 if (1 == domain_bc) {
210 system_l0.
set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
211 system_l0.
set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
212 system_l1.
set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
213 system_l1.
set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
215 system_even.
set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
216 system_even.
set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
217 system_odd.
set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
218 system_odd.
set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
222 system_l0.
set(index01, index01-1) = 1. ;
223 system_l0.
set(index01, index01) = 1. ;
224 system_l1.
set(index01, index01-1) = 1. ;
225 system_l1.
set(index01, index01) = 1. ;
226 system_even.
set(index, index-1) = 1. ;
227 system_even.
set(index, index) = 1. ;
228 system_odd.
set(index, index-1) = 1. ;
229 system_odd.
set(index, index) = 1. ;
231 system_l0.
set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
232 system_l0.
set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
233 system_l1.
set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
234 system_l1.
set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
235 system_even.
set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
236 system_even.
set(index+1, index) = (nr-1)*(nr-1)/alpha ;
237 system_odd.
set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
238 system_odd.
set(index+1, index) = (nr-1)*(nr-1)/alpha ;
246 if ((1 == domain_bc)&&(bc_ced))
247 alpha = -0.25/alpha ;
248 if (1 == domain_bc) {
249 system_l0.
set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
250 system_l1.
set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
252 system_even.
set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
253 system_odd.
set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
257 system_l0.
set(index01, index01) = 1. ;
258 system_l1.
set(index01, index01) = 1. ;
259 system_even.
set(index, index) = 1. ;
260 system_odd.
set(index, index) = 1. ;
262 system_l0.
set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
263 system_l1.
set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
264 system_even.
set(index+1, index) = (nr-1)*(nr-1)/alpha ;
265 system_odd.
set(index+1, index) = (nr-1)*(nr-1)/alpha ;
270 for (
int lz=2; lz<=domain_bc; lz++) {
273 if ((lz == domain_bc)&&(bc_ced))
274 alpha = -0.25/alpha ;
275 system_l0.
set(index01, index01+1) = 1. ;
276 system_l0.
set(index01, index01+2) = -1. ;
277 system_l1.
set(index01, index01+1) = 1. ;
278 system_l1.
set(index01, index01+2) = -1. ;
280 system_even.
set(index, index+1) = 1. ;
281 system_even.
set(index, index+2) = -1. ;
282 system_odd.
set(index, index+1) = 1. ;
283 system_odd.
set(index, index+2) = -1. ;
286 system_l0.
set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
287 system_l0.
set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
288 system_l1.
set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
289 system_l1.
set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
291 system_even.
set(index, index) = -(nr-2)*(nr-2)/alpha ;
292 system_even.
set(index, index+1) = (nr-1)*(nr-1)/alpha ;
293 system_odd.
set(index, index) = -(nr-2)*(nr-2)/alpha ;
294 system_odd.
set(index, index+1) = (nr-1)*(nr-1)/alpha ;
298 if (lz == domain_bc) {
299 system_l0.
set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
300 system_l0.
set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
301 system_l1.
set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
302 system_l1.
set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
304 system_even.
set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
305 system_even.
set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
306 system_odd.
set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
307 system_odd.
set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
311 system_l0.
set(index01, index01-1) = 1. ;
312 system_l0.
set(index01, index01) = 1. ;
313 system_l1.
set(index01, index01-1) = 1. ;
314 system_l1.
set(index01, index01) = 1. ;
315 system_even.
set(index, index-1) = 1. ;
316 system_even.
set(index, index) = 1. ;
317 system_odd.
set(index, index-1) = 1. ;
318 system_odd.
set(index, index) = 1. ;
320 system_l0.
set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
321 system_l0.
set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
322 system_l1.
set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
323 system_l1.
set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
324 system_even.
set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
325 system_even.
set(index+1, index) = (nr-1)*(nr-1)/alpha ;
326 system_odd.
set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
327 system_odd.
set(index+1, index) = (nr-1)*(nr-1)/alpha ;
332 assert(index01 == system01_size) ;
333 assert(index == system_size) ;
338 if (par_mat != 0x0) {
352 const Matrice& sys_even = (need_matrices ? system_even :
354 const Matrice& sys_odd = (need_matrices ? system_odd :
365 int m_q, l_q, base_r ;
366 for (
int k=0; k<np+2; k++) {
367 for (
int j=0; j<nt; j++) {
369 if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q >= l_min)) {
371 int sys_size = (l_q < 2 ? system01_size : system_size) ;
372 int nl = (l_q < 2 ? 1 : 2) ;
377 int nr = mgrid.
get_nr(0) ;
380 if (isexcised==
false){
383 double val_c = 0.; pari = 1 ;
384 for (
int i=0; i<nr-nl; i++) {
385 val_c += pari*coef(0, k, j, i) ;
388 rhs.
set(index) = val_c ;
392 double der_c = 0.; pari = 1 ;
393 for (
int i=0; i<nr-nl-1; i++) {
394 der_c += pari*(2*i+1)*coef(0, k, j, i) ;
397 rhs.
set(index) = der_c / alpha ;
402 for (
int i=0; i<nr-nl; i++) {
403 val_b += coef(0, k, j, i) ;
404 der_b += 4*i*i*coef(0, k, j, i) ;
409 for (
int i=0; i<nr-nl-1; i++) {
410 val_b += coef(0, k, j, i) ;
411 der_b += (2*i+1)*(2*i+1)*coef(0, k, j, i) ;
415 rhs.
set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
419 rhs.
set(index) = -val_b ;
421 rhs.
set(index+1) = -der_b/alpha ;
425 if ((1 == domain_bc)&&(bc_ced))
426 alpha = -0.25/alpha ;
428 int nr_lim = nr - n_conditions ;
429 val_b = 0 ; pari = 1 ;
430 for (
int i=0; i<nr_lim; i++) {
431 val_b += pari*coef(1, k, j, i) ;
434 rhs.
set(index) += val_b ;
437 der_b = 0 ; pari = -1 ;
438 for (
int i=0; i<nr_lim; i++) {
439 der_b += pari*i*i*coef(1, k, j, i) ;
442 rhs.
set(index) += der_b/alpha ;
446 for (
int i=0; i<nr_lim; i++)
447 val_b += coef(1, k, j, i) ;
449 for (
int i=0; i<nr_lim; i++) {
450 der_b += i*i*coef(1, k, j, i) ;
453 rhs.
set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
457 rhs.
set(index) = -val_b ;
458 if (derivative) rhs.
set(index+1) = -der_b / alpha ;
464 if ((1 == domain_bc)&&(bc_ced))
465 alpha = -0.25/alpha ;
467 int nr_lim = nr - 1 ;
469 for (
int i=0; i<nr_lim; i++)
470 val_b += coef(1, k, j, i) ;
472 for (
int i=0; i<nr_lim; i++) {
473 der_b += i*i*coef(1, k, j, i) ;
476 rhs.
set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
480 rhs.
set(index) = -val_b ;
481 if (derivative) rhs.
set(index+1) = -der_b / alpha ;
487 for (
int lz=2; lz<=domain_bc; lz++) {
489 if ((lz == domain_bc)&&(bc_ced))
490 alpha = -0.25/alpha ;
492 int nr_lim = nr - n_conditions ;
493 val_b = 0 ; pari = 1 ;
494 for (
int i=0; i<nr_lim; i++) {
495 val_b += pari*coef(lz, k, j, i) ;
498 rhs.
set(index) += val_b ;
501 der_b = 0 ; pari = -1 ;
502 for (
int i=0; i<nr_lim; i++) {
503 der_b += pari*i*i*coef(lz, k, j, i) ;
506 rhs.
set(index) += der_b/alpha ;
510 for (
int i=0; i<nr_lim; i++)
511 val_b += coef(lz, k, j, i) ;
513 for (
int i=0; i<nr_lim; i++) {
514 der_b += i*i*coef(lz, k, j, i) ;
517 rhs.
set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
521 rhs.
set(index) = -val_b ;
522 if (derivative) rhs.
set(index+1) = -der_b / alpha ;
535 solut = (l_q%2 == 0 ? sys_even.
inverse(rhs) :
539 int dec = (base_r ==
R_CHEBP ? 0 : 1) ;
541 if (isexcised==
false){
543 coef.
set(0, k, j, nr-2-dec) = solut(index) ;
546 coef.
set(0, k, j, nr-1-dec) = solut(index) ;
549 coef.
set(0, k, j, nr-1) = 0 ;
553 for (nl=1; nl<=n_conditions; nl++) {
554 int ii = n_conditions - nl + 1 ;
555 coef.
set(1, k, j, nr-ii) = solut(index) ;
561 coef.
set(1,k,j,nr-1)=solut(index);
564 for (
int lz=2; lz<=domain_bc; lz++) {
566 for (nl=1; nl<=n_conditions; nl++) {
567 int ii = n_conditions - nl + 1 ;
568 coef.
set(lz, k, j, nr-ii) = solut(index) ;
574 for (
int lz=0; lz<=domain_bc; lz++)
575 for (
int i=0; i<mgrid.
get_nr(lz); i++)
576 coef.
set(lz, k, j, i) = 0 ;