First Commit

Etensor/Examples/demo_Etensor.m

+close all
+clear all
+
+%% création des indices: un nom et un range
+i = Eindex('i',2);
+j = Eindex('j',2);
+k = Eindex('k',2);
+l = Eindex('l',2);
+
+%% création de tenseurs
+%ordre 0
+A0 = Etensor(3)
+%%ordre 1 : données (vecteur colonne) et un indice, e.g. i.
+A1 = Etensor(rand(2,1),i)
+B1 = Etensor(rand(2,1),i)
+
+%% changement d'indice
+C1 = B1.idx(j);
+
+%% ordre 2 : une matrice et deux indices
+A2 = Etensor(rand(2,2),i,j)
+B2 = Etensor(rand(2,2),i,j)
+
+%% changement d'indice
+C2 = B2.idx(k,l)
+%% permutation des indices
+D2 = B2.permute(j,i)
+
+%% l'identité
+I = Etensor(eye(2),i,j)
+
+%% quelques opérations
+%A0 * B1_i
+contraction(A0,B1)
+
+%% A1_i * B1_i
+contraction(A1,B1) %ici A1 et B1 dépendent déjà de i
+
+%% A1_i * B1_j
+contraction(A1,B1.idx(j)) %on change d'indice à la volée
+
+%% autre exemple
+contraction(A2.idx(i,j),B2.idx(j,k))
+
+%% autre façon
+contraction(A2.idx(i,j),B2.idx(k,l),I.idx(j,k))
+
+%% en passant par un tenseur d'ordre 4
+tmp = contraction(A2.idx(i,j),B2.idx(k,l))
+contraction(tmp,I.idx(j,k))
+
+%% si on a des objets pour lesquels l'opérateur de différentiation est
+%défini...
+
+A0 = Etensor(sym('x*y+sin(x)'))
+coords = Etensor([sym('x');sym('y')],i)
+A1 = Etensor([sym('x*y+sin(x)');sym('x+y')],i)
+
+%% gradient d'un scalaire
+nabla(A0,coords)
+
+%% gradient d'un vecteur
+nabla(A1,coords.idx(j))
+
+%% divergence
+nabla(A1.idx(l),coords.idx(l))
+
+%% plus complexe
+tmp = contraction(A0,A1)
+tmp = nabla(tmp.idx(i),coords.idx(j))
+contraction(tmp,A1.idx(i))
+
+
+

Etensor/Src/Eindex.m

+%Index class Tensor expression with Einstein notation for the indices
+%see examples for usage.
+%
+% This file is subject to the terms and conditions defined in
+% file 'LICENSE.txt', which is part of this source code package.
+%
+% Principal developer : Adrien Leygue (Adrien.Leygue@ec-nantes.fr)
+%
+
+classdef Eindex < handle
+    properties
+        symbol
+        range
+        value
+    end
+    methods
+        function obj = Eindex(symbol,range)
+            narginchk(2,2);
+            validateattributes(symbol,{'char'},{},mfilename,'symbol',1);
+            validateattributes(range,{'numeric'},{'scalar','integer','finite','nonnegative'},mfilename,'range',2);
+            obj.symbol = symbol;
+            obj.range = range;
+            obj.value = 1;
+        end
+        function disp(obj)
+            a = '';
+            for i=1:numel(obj)
+                a = sprintf('%s %s[1:%d]',a,obj(i).symbol,obj(i).range);
+            end
+            disp(a);
+        end
+        function reset(obj)
+            [obj.value] = deal(1);
+        end
+        function result = getRelevant(obj)
+            result = (obj([obj.range]>1));
+        end
+        function result = ismax(obj)
+            result = all([obj.value]==[obj.range]);
+        end
+        function increment(obj)
+            assert(~ismax(obj),'cannot increment')
+            if numel(obj)==1
+                obj.value = obj.value+1;
+                return;
+            end
+            sz = [obj.range];
+            tmp = cell(1,numel(obj));
+            [tmp{:}] = ind2sub(sz,sub2ind(sz,obj.value)+1);
+            [obj.value] = deal(tmp{:});
+        end
+        function [result,factor] = voigtIndex(obj,sample1)
+            %obj should be of length 2 or 4
+            %range should be 3
+            switch numel(obj)
+                case 4
+                    [result1,factor1] = voigtIndex(obj(1:2),sample1);
+                    [result2,factor2] = voigtIndex(obj(3:4),sample1);
+                    result = [result1(1) result2(1)];
+                    factor = factor1*factor2;
+                case 2
+                    tmp = [obj.value];
+                    if tmp(1)==tmp(2)
+                        result = [tmp(1) 1];
+                        factor = sample1;
+                    elseif sum(tmp)==5
+                        result= [4 1];
+                        factor = sqrt(2*sample1);
+                    elseif sum(tmp)==4
+                        result= [5 1];
+                        factor = sqrt(2*sample1);
+                    elseif sum(tmp)==3
+                        result = [6 1];
+                        factor = sqrt(2*sample1);
+                    else
+                        error('surprising error');
+                    end
+                otherwise
+                    error('voigtIndex only accepts a multiindex of length 2 or 4');
+            end  
+        end
+        
+    end
+end
\ No newline at end of file

FEMTools/Src/FEM_ASSEMBLE.m

+%FEM_ASSEMBLE
+%
+% This file is subject to the terms and conditions defined in
+% file 'LICENSE.txt', which is part of this source code package.
+%
+% Principal developer : Adrien Leygue (Adrien.Leygue@ec-nantes.fr)
+%
+function varargout = FEM_ASSEMBLE(DOM_NAME,COORDS,U_FIELDS,T_FIELDS,K_FIELDS,INTEG_PARAM,varargin)
+    WEAK = [];
+    for i=1:numel(varargin)
+        if ~ischar(varargin{i}), continue; end;
+        switch upper(varargin{i})
+            case 'LHS'
+                WEAK = [WEAK varargin{i+1}]; %#ok<AGROW>
+                continue;
+            case 'RHS'
+                WEAK = [WEAK varargin{i+1}]; %#ok<AGROW>
+                continue;
+            otherwise
+                error(['Unrecognized option:' varargin{i}]);
+        end
+    end
+    
+    %basic check to speed up the assembly
+    sym_galerkin = all(U_FIELDS==T_FIELDS);
+    
+    %gather all fields and check if different fields (handles) have
+    %different names
+    ALL_FIELDS = unique([COORDS; T_FIELDS(:);U_FIELDS(:);K_FIELDS(:)]);
+    assert(numel(unique({ALL_FIELDS.name}))==numel(ALL_FIELDS),'FEM_ASSEMBLE: Different fields share the same name');
+    
+    %Prepare the integration over all the fields
+    Nelem = ALL_FIELDS.set_current_domain(DOM_NAME);
+    INTEG_WEIGHTS = ALL_FIELDS.prepare_integration(INTEG_PARAM);
+    
+    %compute the size of the local and global matrices
+    %compute the assembly indexes in the local matrix
+    [Nglob_col,Nloc_col] = U_FIELDS.compute_offsets;
+    loc_cols = cell(1,numel(U_FIELDS));
+    for i=1:numel(U_FIELDS)
+        loc_cols{i} = U_FIELDS(i).get_local_index;
+    end
+    
+    %sym_galerkin could be used here,... but it would not help much
+    [Nglob_row,Nloc_row] = T_FIELDS.compute_offsets;
+    loc_lines = cell(1,numel(T_FIELDS));
+    for i=1:numel(T_FIELDS)
+        loc_lines{i} = T_FIELDS(i).get_local_index;
+    end
+    
+    %process the weak form to:
+    % 1-extract LHS and RHS
+    % 2-extract a unique element list to compute shape function gradients
+    % 3-extract a catalog of elementary expressions that will be evaluated only once
+    [LHS,RHS,ALL_ELEMENTS,catalog_phi,catalog_value] = process_weak_form(WEAK,U_FIELDS,T_FIELDS,ALL_FIELDS,loc_lines,loc_cols,Nloc_row,Nloc_col);
+    
+    %memory allocation for the assembly step
+    has_rhs = ~isempty(RHS);
+    has_lhs = ~isempty(LHS);
+    if has_lhs
+        sz_loc = Nloc_row*Nloc_col;
+        LHS_i = zeros(Nelem*sz_loc,1);
+        LHS_j = zeros(Nelem*sz_loc,1);
+        LHS_values = zeros(Nelem*sz_loc,1);
+    end
+    if has_rhs
+        RHS_i = zeros(Nelem*Nloc_row,1);
+        RHS_values = zeros(Nelem*Nloc_row,1);
+    end
+    
+    N_cat_phi = size(catalog_phi,1);
+    N_cat_value = size(catalog_value,1);
+    
+    eval_catalog_phi = cell(1,N_cat_phi);
+    eval_catalog_value = zeros(N_cat_value,numel(INTEG_WEIGHTS));
+    
+    rep_row = ones(1, Nloc_row);
+    rep_col = ones(1, Nloc_col);
+    
+    %Loop over the elements
+    for el=1:Nelem
+        
+        %local to global assembly map
+        map_i = T_FIELDS.get_local_to_global_map(el);
+        
+        %should check the use of iJ
+        [J,detJ] = COORDS.jacobian_integration(el);
+        WJ = INTEG_WEIGHTS .* detJ';
+        
+        %iJ = cellfun(@inv,J,'uniformoutput',0);
+        
+        for i=1:N_cat_phi
+            eval_catalog_phi{i} = ALL_ELEMENTS(catalog_phi(i,1)).gen_eval_phi(J,catalog_phi(i,2:end));
+        end
+        for i=1:N_cat_value
+            eval_catalog_value(i,:) = ALL_FIELDS(catalog_value(i,1)).eval_field_integration(el,J,catalog_value(i,4:end),catalog_value(i,2));
+        end
+        
+        if has_lhs
+            A_loc = zeros(Nloc_row*Nloc_col,1);
+            
+            %local to global assembly map
+            if sym_galerkin
+                map_j = map_i;
+            else
+                map_j = U_FIELDS.get_local_to_global_map(el);
+            end
+            
+            for tt = 1:numel(LHS)
+                
+                D = WJ*LHS(tt).factor; %scalar value
+                %aux_fields
+                if LHS(tt).has_known
+                    D = D.* prod(eval_catalog_value(LHS(tt).known_entries,:),1);
+                end
+                
+                tmp = eval_catalog_phi{LHS(tt).test_entry}*diag(D)*eval_catalog_phi{LHS(tt).unknown_entry}';
+                A_loc(LHS(tt).loc_idx)=A_loc(LHS(tt).loc_idx) + tmp(:);
+                
+                
+            end
+            
+            tmp = (1:sz_loc)+(el-1)*sz_loc;
+            LHS_i(tmp) = map_i(rep_col,:)';
+            LHS_j(tmp) = map_j(rep_row,:);
+            LHS_values(tmp) = A_loc(:);
+        end
+        
+        if has_rhs
+            B_loc = zeros(Nloc_row,1);
+            for tt = 1:numel(RHS)
+                D = WJ*RHS(tt).factor; %scalar value
+                %aux_fields
+                if RHS(tt).has_known
+                    D = D.* prod(eval_catalog_value(RHS(tt).known_entries,:),1);
+                end
+                B_loc(RHS(tt).loc_lines)=B_loc(RHS(tt).loc_lines) + eval_catalog_phi{RHS(tt).test_entry}*D';            
+            end
+            tmp = (1:Nloc_row)+(el-1)*Nloc_row;
+            RHS_i(tmp) = map_i(:);
+            RHS_values(tmp) = B_loc(:);
+        end
+        
+    end
+    
+    varargout = cell(1,2);
+    varargout{1} = sparse(Nglob_row,Nglob_col);
+    varargout{2} = sparse(Nglob_row,1);
+    if has_lhs
+        varargout{1} = sparse(LHS_i,LHS_j,LHS_values,Nglob_row,Nglob_col);
+    end
+    if has_rhs
+        varargout{2} = accumarray(RHS_i,RHS_values,[Nglob_row 1]);
+    end 
+    if nargout==1
+        varargout = varargout([has_lhs has_rhs]);
+    end
+    
+end
+

FEMTools/Src/FEM_ASSEMBLE_PC.m

+function varargout = FEM_ASSEMBLE_PC(DOM_NAME,COORDS,U_FIELDS,T_FIELDS,K_FIELDS,INTEG_PARAM,varargin)
+%
+% This file is subject to the terms and conditions defined in
+% file 'LICENSE.txt', which is part of this source code package.
+%
+% Principal developer : Adrien Leygue (Adrien.Leygue@ec-nantes.fr)
+%
+LHS = [];
+RHS = [];
+if ~isempty(K_FIELDS);
+INTEG_VALUES = cell(numel(K_FIELDS),max([K_FIELDS.Ncomp])+1);
+end
+
+for i=1:numel(varargin)
+    if ~ischar(varargin{i}), continue; end;
+    switch upper(varargin{i})
+        case 'LHS'
+            LHS = [LHS varargin{i+1}]; %#ok<AGROW>
+            continue;
+        case 'RHS'
+            RHS = [RHS varargin{i+1}]; %#ok<AGROW>
+            continue;
+        case 'INTEG_VALUES'
+            assert(~isempty(K_FIELDS),'INTEG_VALUES needs known fields');
+            [found,pos] = ismember(varargin{i+1}{1},{K_FIELDS.name});
+            assert(found,'provided fields not found among the known fields');
+            assert(isempty(INTEG_VALUES{pos,1}),'duplicate provided field');
+            %INTEG_VALUES{pos,numel(varargin{i+1})} = [];
+            %INTEG_VALUES(pos,:) = varargin{i+1}; 
+            INTEG_VALUES(pos,1:(K_FIELDS(pos).Ncomp+1)) = varargin{i+1};
+            continue;
+        otherwise
+            error(['Unrecognized option:' varargin{i}]);
+    end
+end
+
+
+%gather all fields and check if different fields (handles) have
+%different names
+ALL_FIELDS = unique([COORDS; T_FIELDS(:);U_FIELDS(:);K_FIELDS(:)]);
+assert(all([ALL_FIELDS.type]~=2),'Value fields not supported');
+assert(numel(unique({ALL_FIELDS.name}))==numel(ALL_FIELDS),'FEM_ASSEMBLE: Different fields share the same name');
+ALL_FIELDS.set_current_domain(DOM_NAME);
+COORDS.precompute_interpolants(DOM_NAME,INTEG_PARAM,ALL_FIELDS);
+
+T_offsets = cumsum([T_FIELDS.Nnodes].*[T_FIELDS.Ncomp]);
+T_offsets = T_offsets(:);
+Nrows = T_offsets(end);
+T_offsets = [0;T_offsets(1:end-1)];
+T_comp_offsets = [T_FIELDS.Nnodes];
+
+if ~isempty(U_FIELDS)
+    U_offsets = cumsum([U_FIELDS.Nnodes].*[U_FIELDS.Ncomp]);
+    U_offsets = U_offsets(:);
+    Ncols = U_offsets(end);
+    U_offsets = [0;U_offsets(1:end-1)];
+    U_comp_offsets = [U_FIELDS.Nnodes];
+end
+
+weights = COORDS.CURRENT_INTERP.PCIntegrationWeights{(COORDS.CURRENT_INTERP.PCIntegrationParam==INTEG_PARAM)};
+
+has_rhs = ~isempty(RHS);
+has_lhs = ~isempty(LHS);
+if has_lhs
+    LHS = LHS(find([LHS.factor]));
+    all_i = [];
+    all_j = [];
+    all_v = [];
+    LHS = process_weak(LHS,U_FIELDS,T_FIELDS,K_FIELDS);
+    for t = 1:numel(LHS)
+        term = LHS(t);
+        tmp = term.factor*weights;
+        U_INTERP = U_FIELDS(term.unknown_id).CURRENT_INTERP;
+        T_INTERP = T_FIELDS(term.test_id).CURRENT_INTERP;
+        for k = 1:numel(term.known_id)
+            id = term.known_id(k);
+            if (~isempty(INTEG_VALUES{id,term.known_comp(k)})) && (term.known_diff_symbol(k)==1)
+                tmp = tmp.* INTEG_VALUES{id,term.known_comp(k)+1}(:);
+            else
+                K_INTERP = K_FIELDS(id).CURRENT_INTERP;
+                tmp = tmp.* (K_INTERP.PCMatrix{find(INTEG_PARAM==K_INTERP.PCIntegrationParam),term.known_diff_symbol(k)}*K_FIELDS(id).get_values('ALL',term.known_comp(k))); %#ok<*FNDSB>
+            end
+        end
+        %tmpA = bsxfun(@times,T_INTERP.PCMatrix{find(INTEG_PARAM==T_INTERP.PCIntegrationParam),term.test_diff_symbol},tmp)'*U_INTERP.PCMatrix{find(INTEG_PARAM==U_INTERP.PCIntegrationParam),term.unknown_diff_symbol};
+        tmpA = T_INTERP.PCMatrix{find(INTEG_PARAM==T_INTERP.PCIntegrationParam),term.test_diff_symbol}'*spdiags(tmp,0,numel(tmp),numel(tmp))*U_INTERP.PCMatrix{find(INTEG_PARAM==U_INTERP.PCIntegrationParam),term.unknown_diff_symbol};
+
+        [tmp_i,tmp_j,tmp_v] = find(tmpA);
+        tmp_i = T_offsets(term.test_id)+ T_comp_offsets(term.test_id)*(term.test_comp-1)+tmp_i;
+        tmp_j = U_offsets(term.unknown_id)+ U_comp_offsets(term.unknown_id)*(term.unknown_comp-1)+tmp_j;
+        
+        all_i = [all_i;tmp_i]; %#ok<AGROW>
+        all_j = [all_j;tmp_j]; %#ok<AGROW>
+        all_v = [all_v;tmp_v]; %#ok<AGROW>
+    end
+    A = sparse(all_i,all_j,all_v,Nrows,Ncols);
+end
+if has_rhs
+    RHS = RHS(find([RHS.factor]));
+    b = zeros(Nrows,1);
+    RHS = process_weak(RHS,U_FIELDS,T_FIELDS,K_FIELDS);
+    for t = 1:numel(RHS)
+        term = RHS(t);
+        T_INTERP = T_FIELDS(term.test_id).CURRENT_INTERP;
+        tmp = term.factor*weights;
+        for k = 1:numel(term.known_id)
+            id = term.known_id(k);
+            if (~isempty(INTEG_VALUES{id,term.known_comp(k)})) && (term.known_diff_symbol(k)==1)
+                tmp = tmp.* INTEG_VALUES{id,term.known_comp(k)+1}(:);
+            else
+                
+                K_INTERP = K_FIELDS(id).CURRENT_INTERP;
+                tmp = tmp.* (K_INTERP.PCMatrix{find(INTEG_PARAM==K_INTERP.PCIntegrationParam),term.known_diff_symbol(k)}*K_FIELDS(id).get_values('ALL',term.known_comp(k)));
+            end
+        end
+        tmpb = tmp'* T_INTERP.PCMatrix{find(INTEG_PARAM==T_INTERP.PCIntegrationParam),term.test_diff_symbol};
+        b_offset = T_offsets(term.test_id)+ T_comp_offsets(term.test_id)*(term.test_comp-1);
+        b((1:numel(tmpb))+b_offset) = b((1:numel(tmpb))+b_offset) + tmpb';
+    end
+end
+
+varargout = cell(1,2);
+varargout{1} = sparse(Nrows,Ncols);
+varargout{2} = sparse(Nrows,1);
+if has_lhs
+    varargout{1} = A;
+end
+if has_rhs
+    varargout{2} = b;
+end
+if nargout==1
+    varargout = varargout([has_lhs has_rhs]);
+end
+end
+
+function weak = process_weak(weak,U_FIELDS,T_FIELDS,K_FIELDS)
+for t = 1:numel(weak)
+    tmp = strcmp(weak(t).test_id{1},{T_FIELDS.name});
+    msg = ['could not find ' weak(t).test_id{1} ' among the provided test fields'];
+    assert(nnz(tmp)==1,msg);
+    weak(t).test_id = find(tmp);
+    weak(t).test_diff_symbol = ELEMENT.diff_id(weak(t).test_diff_symbol);
+    
+    if ~isempty(weak(t).unknown_id)
+        tmp = strcmp(weak(t).unknown_id{1},{U_FIELDS.name});
+        msg = ['could not find ' weak(t).unknown_id{1} ' among the provided unknown fields'];
+        assert(nnz(tmp)==1,msg);
+        weak(t).unknown_id = find(tmp);
+        weak(t).unknown_diff_symbol = ELEMENT.diff_id(weak(t).unknown_diff_symbol);
+    end
+    
+    kn_id = zeros(1,numel(weak(t).known_id));
+    kn_diff = zeros(1,numel(weak(t).known_id));
+    for k = 1:numel(weak(t).known_id)
+        tmp = strcmp(weak(t).known_id{k},{K_FIELDS.name});
+        msg = ['could not find ' weak(t).known_id{k} ' among the provided known fields'];
+        assert(nnz(tmp)==1,msg);
+        kn_id(k) = find(tmp);
+        kn_diff(k) = ELEMENT.diff_id(weak(t).known_diff_symbol(:,k));
+    end
+    weak(t).known_id = kn_id;
+    weak(t).known_diff_symbol = kn_diff;
+end
+end
+

FEMTools/Src/FEM_ASSEMBLE_PGD.m

+%FEM_ASSEMBLE_PGD
+%
+% This file is subject to the terms and conditions defined in
+% file 'LICENSE.txt', which is part of this source code package.
+%
+% Principal developer : Adrien Leygue (Adrien.Leygue@ec-nantes.fr)
+%
+function varargout = FEM_ASSEMBLE_PGD_2(DOM_NAME,COORDS,U_FIELDS,T_FIELDS,K_FIELDS,INTEG_PARAM,varargin)
+    
+    old_school = false;
+    WEAK = [];
+    for i=1:numel(varargin)
+        if ~ischar(varargin{i}), continue; end;
+        switch upper(varargin{i})
+            case 'LHS'
+                WEAK = [WEAK varargin{i+1}]; %#ok<AGROW>
+                continue;
+            case 'RHS'
+                WEAK = [WEAK varargin{i+1}]; %#ok<AGROW>
+                continue;
+            case 'OLD_SCHOOL'
+                old_school = true;
+                continue;
+            otherwise
+                error(['Unrecognized option:' varargin{i}]);
+        end
+    end
+    %basic check to speed up the assembly
+    %sym_galerkin = all(U_FIELDS==T_FIELDS);
+    
+    %gather all fields and check if different fields (handles) have
+    %different names
+    ALL_FIELDS = unique([COORDS; T_FIELDS(:);U_FIELDS(:);K_FIELDS(:)]);
+    assert(numel(unique({ALL_FIELDS.name}))==numel(ALL_FIELDS),'FEM_ASSEMBLE: Different fields share the same name');
+    
+    %Prepare the integration over all the fields
+    Nelem = ALL_FIELDS.set_current_domain(DOM_NAME);
+    INTEG_WEIGHTS = ALL_FIELDS.prepare_integration(INTEG_PARAM);
+    
+    %compute the size of the local and global matrices
+    %compute the assembly indexes in the local matrix
+    [Nglob_col,Nloc_col] = U_FIELDS.compute_offsets;
+    loc_cols = cell(1,numel(U_FIELDS));
+    for i=1:numel(U_FIELDS)
+        loc_cols{i} = U_FIELDS(i).get_local_index;
+    end
+    
+    %sym_galerkin could be used here,... but it would not help much
+    [Nglob_row,Nloc_row] = T_FIELDS.compute_offsets;
+    loc_lines = cell(1,numel(T_FIELDS));
+    for i=1:numel(T_FIELDS)
+        loc_lines{i} = T_FIELDS(i).get_local_index;
+    end
+    
+    %process the weak form to:
+    % 1-extract LHS and RHS
+    % 2-extract a unique element list to compute shape function gradients
+    % 3-extract a catalog of elementary expressions that will be evaluated only once
+    
+    %TODO : ajouter les facteurs et faire l'expansion des aux-fields
+    [LHS,RHS,ALL_ELEMENTS,catalog_phi,catalog_value] = process_weak_form(WEAK,U_FIELDS,T_FIELDS,ALL_FIELDS,loc_lines,loc_cols,Nloc_row,Nloc_col);
+    has_rhs = ~isempty(RHS);
+    has_lhs = ~isempty(LHS); 
+    
+    %virtual assembly and memory allocation
+    if has_lhs
+        size_term = zeros(1,numel(LHS));
+        for tt = 1:numel(LHS)
+            myterm = LHS(tt);
+            my_i = myterm.test_entry; %test field number
+            my_j = myterm.unknown_entry; %unnown field number
+            size_term(tt) = ALL_ELEMENTS(catalog_phi(my_i,1)).Nphi * ALL_ELEMENTS(catalog_phi(my_j,1)).Nphi;
+        end
+        
+        %memory allocation for virtual assembly
+        LHS_i = zeros(Nelem*sum(size_term),1);
+        LHS_j = zeros(Nelem*sum(size_term),1);
+       
+        idx_lhs = 0;
+        for el = 1:Nelem
+            for tt = 1:numel(LHS)
+                myterm = LHS(tt);
+                    
+                map_i = U_FIELDS(myterm.test_pos).get_detailled_map(el,myterm.test_comp);
+                map_j = T_FIELDS(myterm.unknown_pos).get_detailled_map(el,myterm.unknown_comp);
+                
+                tmp_i = zeros(numel(map_i),numel(map_j));
+                tmp_j = zeros(numel(map_i),numel(map_j));
+                
+                for k = 1:numel(map_i)
+                    tmp_j(k,:) = map_j;
+                end
+                for k = 1:numel(map_j)
+                    tmp_i(:,k) = map_i;
+                end
+                tmp = idx_lhs+(1:(numel(map_i)*numel(map_j)));
+                
+                LHS_i(tmp) = map_i(ones(1,numel(map_j)),:)';
+                LHS_j(tmp) = map_j(ones(1,numel(map_i)),:);
+                idx_lhs = idx_lhs + numel(tmp);
+            end
+        end
+        A = sparse(LHS_i,LHS_j,ones(size(LHS_i)),Nglob_row,Nglob_col);
+        [pattern_i,pattern_j] = find(A);
+        A_pattern = sparse(pattern_i,pattern_j,1:numel(pattern_i),Nglob_row,Nglob_col);
+    end
+    
+    [LHS,RHS,catalog_value] = expand_separated_terms([LHS(:);RHS(:)],catalog_value,ALL_FIELDS);
+    
+    if has_lhs
+        size_term = zeros(1,numel(LHS));
+        for tt = 1:numel(LHS)
+            myterm = LHS(tt);
+            my_i = myterm.test_entry; %test field number
+            my_j = myterm.unknown_entry; %unnown field number
+            size_term(tt) = ALL_ELEMENTS(catalog_phi(my_i,1)).Nphi * ALL_ELEMENTS(catalog_phi(my_j,1)).Nphi;
+        end 
+        LHS_i = zeros(Nelem*sum(size_term),1);
+        LHS_j = zeros(Nelem*sum(size_term),1);
+        LHS_values = zeros(Nelem*sum(size_term),1);
+    end
+    
+    if has_rhs
+        size_term = zeros(1,numel(RHS));
+        for tt = 1:numel(RHS)
+            myterm = RHS(tt);
+            my_i = myterm.test_pos;
+            size_term(tt) = numel(myterm.loc_lines);%ALL_ELEMENTS(my_i).Nphi;
+        end
+        RHS_i = zeros(Nelem*sum(size_term),1);
+        RHS_j = zeros(Nelem*sum(size_term),1);
+        RHS_values = zeros(Nelem*sum(size_term),1);
+    end
+
+    N_cat_phi = size(catalog_phi,1);
+    N_cat_value = size(catalog_value,1);
+    
+    eval_catalog_phi = cell(1,N_cat_phi);
+    eval_catalog_value = zeros(N_cat_value,numel(INTEG_WEIGHTS));
+    
+    %Loop over the elements