Coefficient that varies on domain

General Discussion
1

Hello, I’m trying to implement a coefficient value that varies depending on the area of the domain, this is the definition

//parameters
func real kappa() {
 int ni = 9, nj = 9;
 real r1=1./32.*3./ni, k1=1., k3=1.;
 real r2=1./16.*3./ni, k2=10000.0;

 for(int j=0; j<(nj-1); j++){
  real yj=(2.*j+1.)/(2.*nj);
  for(int i=0; i<(ni-1); i++){
   real xi=(2.*i+1.)/(2.*ni);
   if( (x>=(xi-r1)) && (x<=(xi+r1)) && (y>=(yj-r1)) && (y<=(yj+r1)))
    return k1;
   if ( (x>=(xi-r2)) && (x<=(xi+r2)) && (y>=(yj-r2)) && (y<=(yj+r2)))
    return k2;
  }
 }
return k3;

}

I need this coefficient for the integral for the Darcy Model

   varf ahK(etaK, varphi)
    =  int2d(Th[K])(kappa*grad(etaK)'*grad(varphi));

but I’m getting the following error

varf ahK(etaK, varphi)
  289 :     =  int2d(Th[K])(kappa*grad(etaK)   [dx(etaK), dy(etaK)]'* error operator *  <11Polymorphic>, <10LinearCombI7MGauche4C_F0E>

can someone help me?

2

this is an example of the code I want to implement on FreeFem++

local function func_kappa(x,y)
	ni=27; nj=27;
	-- (not so) aligned to the mesh
	r1=1./32.*3./ni ; k1=1.; k3=1.; -- radius and kappa of the inner squared cell
	r2=1./16.*3./ni ; k2=10000.0; -- radius and kappa of the outer square cell

	for j=0,(nj-1) do
		yj=(2.*j+1.)/(2.*nj);
		for i=0,(ni-1) do
			xi=(2.*i+1.)/(2.*ni);
			if ( (x>=(xi-r1)) and (x<=(xi+r1)) and (y>=(yj-r1)) and (y<=(yj+r1)))  then
				return k1;
			end
			if ( (x>=(xi-r2)) and (x<=(xi+r2)) and (y>=(yj-r2)) and (y<=(yj+r2)))  then
				return k2;
			end
		end
	end
	return k3; -- kappa of the embedding matrix
end


kappa = func_kappa,
3

Did that work? If not can you upload your entire code?

4

The technics is very slow, you function kappa is defined by
with a Gh grid (x_i,y_i) of size ni, nj you want

let us call dist00(x,y) is the minimal distance of (x,y) the the point of the grid in norme infini
now the fonction is just k3+(dist00(x,y) < r1)(k1-k3) + (dist00(x,y) < r2)(k2-k1-k3)

the trick to code dist00(x,y) arrive to morrow!

because

5

It worked like this:

//parameters
func real varkappa(int in, int jn){
  int ni = in, nj = jn;
  real r1=1./32.*3./ni, k1=1., k3=1.;
  real r2=1./16.*3./ni, k2=10000.0;

  for(int j=0; j<(nj-1); j++){
    real yj=(2.*j+1.)/(2.*nj);
    for(int i=0; i<(ni-1); i++){
      real xi=(2.*i+1.)/(2.*ni);
      if( (x>=(xi-r1)) && (x<=(xi+r1)) && (y>=(yj-r1)) && (y<=(yj+r1))){
        return k1;
      }
      if ( (x>=(xi-r2)) && (x<=(xi+r2)) && (y>=(yj-r2)) && (y<=(yj+r2))){
        return k2;
      }
    }
  }
  return k3;
}

real kappa = varkappa(9,9);
6

so my final version (the complexity of function Kappa is O(1) not O(ni*nj)

yyy.edp (1.1 KB)

1 Like